Bivariate network meta‐analysis for surrogate endpoint evaluation

Surrogate endpoints are very important in regulatory decision making in healthcare, in particular if they can be measured early compared to the long‐term final clinical outcome and act as good predictors of clinical benefit. Bivariate meta‐analysis methods can be used to evaluate surrogate endpoints and to predict the treatment effect on the final outcome from the treatment effect measured on a surrogate endpoint. However, candidate surrogate endpoints are often imperfect, and the level of association between the treatment effects on the surrogate and final outcomes may vary between treatments. This imposes a limitation on methods which do not differentiate between the treatments. We develop bivariate network meta‐analysis (bvNMA) methods, which combine data on treatment effects on the surrogate and final outcomes, from trials investigating multiple treatment contrasts. The bvNMA methods estimate the effects on both outcomes for all treatment contrasts individually in a single analysis. At the same time, they allow us to model the trial‐level surrogacy patterns within each treatment contrast and treatment‐level surrogacy, thus enabling predictions of the treatment effect on the final outcome either for a new study in a new population or for a new treatment. Modelling assumptions about the between‐studies heterogeneity and the network consistency, and their impact on predictions, are investigated using an illustrative example in advanced colorectal cancer and in a simulation study. When the strength of the surrogate relationships varies across treatment contrasts, bvNMA has the advantage of identifying treatment comparisons for which surrogacy holds, thus leading to better predictions.


Introduction
Surrogate endpoints are very important in the drug development process, at both the trial design and the evaluation stage. They are particularly useful when they can provide early measurement of the treatment effect, in settings where a long follow up time is required before measurement of the final clinical outcome [6]. This is often the case in cancer where overall survival (OS) is of primary interest whilst other outcomes such as progression-free survival (PFS) potentially can be used to measure the effect of a treatment earlier. Alternatively, PFS may be of primary interest and tumour response (TR) is then investigated as a short term surrogate endpoint to PFS. Before they can be used in evaluation of new health technologies, candidate surrogate endpoints have to be assessed for their predictive value of the treatment effect on the final clinical outcome. Surrogate outcomes are validated by estimating the pattern of association between the treatment effects on surrogate and final endpoints across trials using meta-analytic techniques [10,8,7,18,3,2].
Multivariate meta-analysis methods are used to obtain average treatment effects on multiple endpoints while taking account of the correlations between them [20,19,22,4] and, as such, are suitable tools for modelling surrogate endpoints [3,2]. Bivariate meta-analysis (bvMA) of treatment effects on a surrogate and a final outcome allows for both the validation of a surrogate endpoint 1 arXiv:1807.08928v1 [stat.ME] 24 Jul 2018 and for making predictions of an unobserved treatment effect on the final clinical outcome from observed treatment effects on a surrogate endpoint.
Candidate surrogate endpoints often are not perfect, and the association patterns between the treatment effects on the surrogate and final outcomes may vary between treatments. Whilst bvMA, described in detail in Section 3, can be used to model surrogate endpoints, it does not differentiate between the treatment options. This is an important issue when the evidence base consists of multiple trials of different treatments in different populations. Network meta-analysis (NMA) combines data from trials investigating heterogeneous treatment contrasts and has the advantage of estimating effects for all treatment contrasts individually. At the same time, the consistency assumption allows for combining evidence (and borrowing strength) across the treatment contrasts. In this paper, we will exploit this framework to model surrogacy relationships by extending previously reported methods of multivariate NMA [1,11,12,13].
In the multivariate NMA (mvNMA), true treatment effects on multiple correlated outcomes are assumed to follow separate multivariate distributions for each treatment contrast. The previously reported mvNMA methods typically simplified the between-studies variance-covariance structures by assuming homogeneity of the correlations and the heterogeneity parameters across the treatment contrasts. We relax this assumption of homogeneity to model in detail different association patterns between the effects on surrogate and final endpoints separately for different treatments or treatment classes. By allowing such patterns to differ across treatments, the methodology can help to identify treatments for which the surrogacy holds and to improve predictions. The network consistency assumption provides a framework for combining evidence across treatment contrasts and hence modelling and distinguishing surrogacy patterns within and across treatment contrasts. These two levels of surrogacy enable predictions of the treatment effect on the final outcome in a new study investigating either an existing treatment in a new population (within treatment contrast surrogacy) or a new treatment (across treatment surrogacy). We extend the second order consistency condition, described by Lu and Ades [16] in a univariate NMA, to the bivariate case, in order to gain additional precision in modelling surrogate relationships.
We illustrate the use of the methods, described in Sections 3 and 4, by fitting them to data from an example in advanced colorectal cancer introduced in Section 2 (results presented in Section 5). To illustrate the motivation and the application of the methodology in a more detailed and controlled manner, we apply the methods to simulated data scenarios (Section 6). We fit all models in a Bayesian framework using the OpenBUGS software.

Illustrative example
We use data from randomised controlled trials (RCTs) in advanced colorectal cancer (aCRC), investigating a range of different treatment options, to illustrate how the hierarchical bvNMA models, proposed in this paper, can differentiate the association patterns between the treatment contrasts or classes whilst borrowing strength across treatment contrasts. Data were obtained from four published systematic reviews of RCTs investigating pharmacological treatments in aCRC, categorised into classes with respect to their mechanism of action. These were targeted therapies including antiangiogenic treatments targeting vascular endothelial growth factor (anti-VEGF) [21], anti epidermal growth factor receptor inhibitors (EGFRi) [9], humanised monoclonal antibody targeting integrin receptors (anti-IgG2) and monoclonal antibody targeting the type 1 insulin-like growth factor receptor (anti-IGF1R) [17] or chemotherapies compared to the targeted therapies [17,14] and combinations of these therapies. 15 RCTs investigated use of anti-VEGF with chemotherapy vs. chemotherapy alone, 24 RCTs of EGFRi with chemotherapy vs. chemotherapy alone, 4 RCTs of EGFRi with chemotherapy vs. anti-VEGF with chemotherapy, 4 RCTs of EGFRi with anti-VEGF and chemotherapy vs. anti-VEGF with chemotherapy, one study (in two subgroups of population) of anti-IgG2 with EGFRi and chemotherapy vs. EGFRi with chemotherapy, one study of anti-IGF1R with chemotherapy vs. chemotherapy alone and one of EGFRi with anti-VEGF and chemotherapy vs. chemotherapy alone. The treatments are summarised in the network diagram of Figure 1. Tumour response (TR) is used as an example of a potential surrogate endpoint to progression free survival (PFS). The models are applied to data representing treatment effect on the two outcomes on log odds ratio (OR) scale for TR and log hazard ratio (HR) scale for PFS. The scatter plot in Figure 1 illustrates the association patterns between the treatment effects on the two outcomes for each treatment contrast. The bivariate random effects meta-analysis (BRMA) model for correlated and normally distributed treatment effects on two outcomes Y 1i and Y 2i is usually presented in the form described by van Houwelingen et al. [20] and Riley et al. [19]: In this model, the treatment effects on the surrogate endpoint Y 1i and on the final outcome Y 2i are assumed to estimate the correlated true effects µ 1i and µ 2i with corresponding within-study variances σ 2 1i and σ 2 2i of the estimates and the within-study correlation ρ wi between them. In this hierarchical framework, these true study-level effects follow a bivariate normal distribution with means (β 1 , β 2 ) corresponding to the two outcomes, between-studies variances τ 2 1 and τ 2 2 and a between-studies correlation ρ. Equation (1) represents the within-study model and (2) is the between-studies model. To implement the model in the Bayesian framework, prior distributions are placed on the mean effects, for example vague prior distributions β 1 ∼ N (0, 10 4 ), β 2 ∼ N (0, 10 4 ) and on the between-studies variances and correlation.
In the general case, for any number of outcomes, a prior distribution has to be placed on the whole variance-covariance matrix or the correlation matrix in such a way to ensure that the variancecovariance matrix is positive semi-definite. This can be achieved by placing an inverse Wishart prior distribution on the variance-covariance matrix [22] or by using a separation strategy, with spherical [16,22] or Cholesky decomposition [22] of the correlation matrix. Alternatively, a product normal formulation of the between-studies model can be used where the model is parameterised in the form of a series of univariate conditional distributions, ensuring that the relationships between the parameters of the model (regression coefficients and conditional variances) and the elements of the between-studies variance-covariance matrix result in the positive semi-definite between-studies variance-covariance matrix [4,3,2].
In the bivariate case, such as considered here, positive semi-definiteness can be achieved by placing prior distributions on the variances, which are restricted to plausible positive values, for example by placing uniform prior distributions on the corresponding standard deviations τ j ∼ U (0, 2), and by choosing a prior distribution for the correlation which restricts it to values between −1 and +1. Here we use a beta distribution to construct such a prior: ρ+1 2 ∼ Beta(1.5, 1.5), as in Burke et al. [5], with details also in the Appendix A.

Bivariate network meta-analysis (bvNMA)
The BRMA model is typically used to obtain average effects on two correlated outcomes, for studies of the same treatment or treatment class. In the context of surrogate endpoints, data on a range of treatments are typically used. BRMA assumes exchangeability of the treatment effects from all studies regardless of the treatment contrast, by assuming that the true effects follow a common (here bivariate normal) distribution (as in Eq. (2)). This model works well for strong surrogate relationships across all treatments. However, the association pattern between treatment effects measured on the surrogate and final outcomes may depend on the treatment contrast. Network meta-analysis can differentiate between treatment contrasts and in bivariate form can be applied to modelling such surrogate relationships within and across treatment contrasts. We first discuss general bvNMA models in Section 4.1 which describe the within-treatment surrogate relationships in the network. In Section 4.2 we extend these models to allow for an additional level of surrogacy across treatment contrasts.

Model 1a: bvNMA
To model correlated treatment effects on a surrogate endpoint and a final outcome for which the correlation varies according to the treatment contrasts, the assumption made by BRMA, that the true effects follow a common distribution, can be replaced by an assumption that the true effects corresponding to different treatment contrasts follow separate distributions. This naturally leads to relaxing the assumption of homogeneity of the between-studies covariance matrix, allowing its elements (the between-studies correlations ρ 1kl,2kl between the treatment effects l vs. k on the surrogate (1) and final (2) outcomes and the heterogeneity parameters for the treatment effects on the two outcomes τ 2 1kl and τ 2 2kl ) to vary across the treatment contrasts kl. To take into account the network structure of the data, we model the treatment effect differences Y jkli between treatments k and l in study i for outcome j = 1, 2, as follows: where k and l denote baseline (control) and experimental treatments respectively in a study i, µ jkli denote the random true treatment effects (differences between the effects of treatments k and l) on outcome j in study i and, and the d jkl are mean treatment effect differences between treatments k and l for each outcome j. We use the first-order consistency assumption, as described by Lu and Ades [16], extended here to the bivariate case. For any three treatments (b, k, l), the treatment differences (µ jkli ) satisfy the following transitivity relations Taking the expectation on both sides gives the consistency equations for the first-order moments which represent the relationships between the treatment contrasts in the population. When b = 1 is a common reference treatment in the network, the treatment effects of each treatment k in the network relative to this common reference treatment 1; the d j,1k are referred to as basic parameters for each outcome j, with d j,11 = 0 and the others are given prior distributions: Prior distributions are also placed on the elements of the between-studies variance-covariance matrix. Similarly as in BRMA, prior distributions for the heterogeneity parameters are selected in a way to ensure that they are restricted to plausible positive values, such as τ jkl ∼ unif (0, 2) and for the correlations to ensure restriction to the values between −1 and 1, e.g. ρ kl +1 2 ∼ Beta(1.5, 1.5), thus guaranteeing that the variance-covariance matrix is positive semi-definite.

Model 1b: bvNMA with second order consistency
The above meta-analytic model 1a assumes consistency of treatment effects on both outcomes. This implies some constraints on the between-studies variance-covariance matrices which can be explicitly introduced to the bvNMA model (3)-(4) by assuming the consistency of the second-order moments. To do so, we extend the approach proposed by Lu and Ades [16] to the bivariate case by taking variances on both sides of the transitivity equation (5), which gives leading to the following relationship between the variances for any three treatments (b, k, l) and for both outcomes j = 1, 2: which gives the second-order consistency conditions (triangle inequalities): In addition, the following condition applies to the covariances: which implies further constraints that are more complex than those in Eq. (10).
To ensure that prior distributions for heterogeneous variance-covariance matrices are appropriate, i.e. to maintain the second-order consistency condition for any three treatments in the network, bivariate ancillary parameters are used, extending the univariate model by To ensure that prior distributions for heterogeneous variance-covariance matrices are appropriate, i.e. to maintain the second-order consistency condition for any three treatments in the network, bivariate ancillary parameters are used, similarly as by Lu and Ades [16], allowing the between-studies variance-covariance matrices to be represented as where γ 2 jk and γ 2 jl are the ancillary parameters: the variances of two correlated random effects ζ jki and ζ jli corresponding to treatment arms k and l (for each outcome j = 1, 2). Prior distributions for the set of between-studies standard deviations τ jkl for each outcome j and each pair of treatments k and l can be given by constructing a prior distribution for a covariance matrix Γ composed of the standard deviations γ jk and the correlations ξ jk,j l between the effects ζ jki and ζ j li , for j, j = 1, 2 and k, l = 1, . . . , n t , where n t is the number of treatments in the network. For the set of values of the elements of matrix Γ to give a resulting set of standard deviations τ jkl that satisfy the second-order consistency rules (10) and (11), the matrix Γ has to be positive semi-definite. This is achieved using a separation strategy with a Cholesky decomposition: Γ = V 1/2 RV 1/2 , where V 1/2 is a 2n t × 2n t diagonal matrix of the standard deviations γ 11 , γ 21 , . . . , γ 1nt , γ 2nt and R is a positive semi-definite 2n t × 2n t matrix of correlations ξ jk,ĵl (block matrix consisting of n t × n t blocks that are of 2 × 2 dimension). Matrix R is represented as R = L T L with L being a 2n t × 2n t upper triangular matrix. To obtain the elements of the matrix L, we extended the method by Wei and Higgins (2013) who describe it for a four dimensional matrix case. Prior distributions are placed on the standard deviations, which need to be restricted to positive values, for example γ j,k ∼ unif (0, 2). The prior distributions placed on the ancillary variances and correlations give implied prior distributions on the between-studies correlations and standard deviations through the formulae (12).

Model 1c: bvNMA with second order consistency and similarity of the ancillary parameters
To ensure second order consistency of the treatment effects in the network, model 1b requires estimation of relatively many parameters. Where data for each treatment k is only available from a limited number of studies, it may be difficult to estimate the individual variances γ 2 1k . To overcome this issue, exchangeability of the ancillary variances may be assumed by replacing individual prior distributions for these parameters with a common distribution; 4.1.4 Model 1d: bvNMA assuming homogeneity of the between-studies variancecovariance matrix When networks are sparse, instead of assuming similarity between the heterogeneity parameters, as in Eq. (13), homogeneity of the between-studies variance-covariance matrix can be assumed, which is a common assumption in the multivariate network meta-analysis or the univariate network meta-analysis of multi-arm trials [15]. The between-studies equation (4) then becomes where µ jkli is the random treatment effect difference between treatments k and l on outcome j in study i and d j,1k is the mean treatment effect difference between treatment k and the reference treatment 1 for outcome j, with d j,11 = 0. A prior distribution is placed on each basic parameter, d j,1k ∼ N (0, 10 3 ). A prior distribution is also placed on the elements of the common betweenstudies variance-covariance matrix: the between-studies standard deviations, τ j ∼ U nif (0, 2) and the correlation, ρ+1 2 ∼ Beta(1.5, 1.5).

Surrogacy criteria:
When using the bivariate NMA models 1a-c to evaluate surrogate endpoints within treatment contrast kl, perfect surrogacy means that This condition is equivalent to assuming that in the linear relationship between true treatment effects on the final and surrogate endpoints, the intercept is zero, µ 2kli = const × µ 1kli and the conditional variance (of the treatment effect on the final outcome conditional on the effect on the surrogate endpoint) is also zero, as in the standard surrogacy models [10,3] (for more details see Appendix B). This criterion describes the between-studies (and populations) surrogacy relationships within the treatment contrasts. The network structure, with the unique surrogate relationship for each treatment contrast, can help us to disentangle information about a surrogate relationship for a particular treatment and to make better predictions, in particular when these relationships are clearly distinct, as illustrated in Section 6. These individual surrogacy relationships enable predicting the treatment effect on the final outcome from the treatment effect measured on a surrogate endpoint in a new study investigating an existing treatment either in a new population or with a different (but already present in the network) control treatment.
When evaluating surrogate endpoints with model 1d, similarly as for models 1a-c the random effects are assumed to follow separate distributions for studies evaluating different treatment contrasts kl but the between-studies correlations and heterogeneity parameters are assumed to be the same across the treatment contrasts. In this case, the perfect surrogacy means that ρ = ±1 and µ 1kli = 0 ⇔ µ 2kli = 0 (16) which differs from the criteria for the models 1a-c with respect to the common correlation across the treatment contrasts. The strength of the association between the treatment effects on the surrogate endpoint and the final clinical outcome will be the same across the treatment contrasts when using model 1d, similarly as when using a pairwise meta-analysis such as BRMA. However, assuming that the random effects for different treatment contrasts follow separate distributions in 6 the bvNMA, we model these association patterns in more detail compared to when using BRMA. In addition, taking into account the network structure of the data results in borrowing of strength across the treatment contrasts when evaluating surrogate relationships.

bvNMA models with borrowing of strength across treatment contrasts
Models 1a-d describe the relationships between the treatment effects on correlated outcomes when for each treatment contrast in the network there is at least one study reporting the treatment effects on both outcomes. In the situation where we want to predict the treatment effect on the final outcome from the effect on the surrogate endpoint in a new study evaluating a new treatment, there will be only this one study reporting the effect for that treatment and only on the surrogate endpoint. In this case the above models 1a-d will result in an estimate of the average effect on the final outcome for the new treatment based solely on the prior distribution, because no data are available for the effect of this treatment on the final outcome. For example, as in the network depicted in Figure 2, if we select treatment A as the reference treatment 1, then the average treatment effect for the new treatment D on the final outcome equals d 2,CD = d 2,AD − d 2,AC and the estimate of the basic parameter d 2,AD would be based on the prior distribution only. As a result, the estimate of the between-studies correlation ρ 1CD,2CD will be based solely on the prior distribution. Therefore the predicted treatment effect on the final outcome from the treatment effect on the surrogate endpoint for the new treatment will not be meaningful. The models 1a-d are designed to evaluate surrogacy patterns across different populations within a treatment contrast and make predictions of a treatment effect on the final outcome in a new study investigating an existing treatment only, but perhaps in a new population.

Model 2a: assuming exchangeability of treatments
To overcome the above issue, an exchangeability assumption can be made to allow for the relationships between the average effects on the two outcomes to be similar across treatment options.
To achieve this for model 1a (Eqs (3), (4) and (6)), instead of placing a prior distribution on each basic parameter (Eq. (7)), it is assumed that the pooled effects for each treatment arm k, θ jk on the two outcomes j = 1, 2, are exchangeable and correlated: for k = 1, . . . , n t . Assuming exchangeability of the effects in each treatment arm θ jk , rather than of the basic parameters d j1k which are the effects relative to a common reference treatment, ensures that the prior distributions for each pair (θ 1k , θ 2k ) are independent. The modelled data Y jkli are the treatment effect differences, thus the ancillary parameters θ jk are not identifiable. However, we are not interested in estimating these parameters, but in estimation of the association between the average treatment effect differences on the two outcomes: d jkl , j = 1, 2. The above formulae imply the association between the average effects: k = 2, . . . , n t , cancelling out the mean values η j . Therefore the mean values can be set arbitrarily to a constant, for example η j = 0. Prior distributions are placed on the elements of the covariance matrix: ω j ∼ U nif (0, 2) and ρt+1 2 ∼ Beta(1.5, 1.5). Borrowing of strength across treatments by relating these parameters enables the model to predict the effect of a new treatment on the final outcome, as in the network scenario in Figure 2.

Model 2b: assuming exchangeability of treatments, with second order consistency constraints
We extend model 1b, described by the general mvNMA (3)-(4) along with the first (6) and second order (10-11) consistency assumptions, to allow for the exchangeability of basic parameters. Similarly as in model 2a, this replaces the individual prior distributions on the basic parameters (7) with the formulae (17)-(18), with prior distributions as in Section 4.2.1.

Model 2c: assuming exchangeability of treatments, with second order consistency constraints and exchangeable variances in treatment arms
This model is an extension of model 1c, assuming similarity across the ancillary parameters described by (13), to allow for the exchangeability of the basic parameters using the formulae (17)- (18), as in models 2a and 2b.

Model 2d: assuming exchangeability of treatments and homogeneity of the between-studies variance-covariance matrix
To extend model 1d to allow for the exchangeability of the correlated basic parameters, we used equations (3) and (14) together with the exchangeability model (17)-(18) on the basic parameters. As in model 1d, prior distributions are placed on the elements of the common between-studies variance-covariance matrix, the between-studies standard deviations, τ j ∼ U nif (0, 2) and the correlation, ρ+1 2 ∼ Beta(1.5, 1.5).

Surrogacy criteria
Similarly as for models 1a-d, the surrogate relationship between the treatment effects on the surrogate endpoint and the final clinical outcome is perfect when the correlation ρ kl = ±1 (ρ ± 1 for model 2d) and no treatment effect on the surrogate endpoint (µ 1kli = 0) will imply no effect on the final outcome (µ 2kli = 0) in the same study i, as described by the formulae (15) and (16). This is the surrogacy relationship between studies (or populations) described within each treatment contrast kl and it is assumed to differ across treatment contrasts. In addition to this, another level of surrogate relationship is described by models 2a-d, the across-treatments surrogacy. Such a surrogate relationship would be perfect if a zero average treatment effect difference between treatments k and l on the surrogate endpoint will imply zero average treatment effect also on the final clinical endpoint for the same treatment contrast. This surrogacy relationship allows prediction of the treatment effect on the final outcome from the treatment effect measured on a surrogate endpoint in a new study investigating a new treatment.

Summary and discussion of models
Bivariate models for network meta-analysis described in this Section differ in their specific assumptions about the heterogeneity between studies and treatments. All components of the between-study models for all models, 1a-d and 2a-d, are summarised in Table 1. The models allow the correlations between the treatment effects on the surrogate and final endpoints to vary to a different degree, and which model is applied in practice can be determined based on the available evidence. The models reduce to the standard meta-analysis model for surrogate endpoints, such as the BRMA model (1)-(2), in a special case of data structure when there are only two treatments in the network, as detailed in Appendix C.

Strategies for model comparison
To make comparisons between the models, in the first instance the between-studies correlations (as well as the mean effects and the heterogeneity parameters) are obtained and compared across the models. Following this, predicted values (and corresponding credible intervals (CrIs)) of the treatment effects on the final outcome are compared to the observed estimates (and corresponding confidence intervals (CIs)) in take-one-out cross-validation procedure. In one study at a time, the estimate of the treatment effect on the final outcome Y 2i is removed (and treated as missing at random) and then this treatment effect is predicted from the treatment effect on the surrogate endpoint, conditional of the data on both outcomes from all the remaining studies in the meta-analysis. The standard deviation of the predicted effectŶ 2i is equal to σ 2 2i + var(μ 2(kl)i |Y 1i , σ 1i , Y 1(−i) , Y 2(−i) ), where Y 1(−i) and Y 2(−i) denote the data from the remaining studies without the validation study i (the treatment contrast subscripts, present only in the network meta-analysis, have been dropped here).
The models were compared with respect to the predicted effects by investigating a number of statistics obtained from each model: a) proportion of the CI of the observed estimate of the effect on the final outcome overlapping with the CrI of the predicted effect, p overlap , b) mean absolute difference between observed estimate of mean effect and predicted mean effect on the final outcome across studies, |m obs − m pred |, c) ratio of the width of intervals; the width of CrI of the predicted effect to the width of the CI of the observed effect (the estimate), w pred /w obs , d) percentage reduction in uncertainty measured by the width of the CrI of the predicted effect when using a mvNMA model, w N M A , compared to the width of the CrI obtained from BRMA w BRM A ; , and e) a new statistic measuring overall performance of a model giving a higher score for models resulting in large p overlap with penalty for overly inflated predicted interval: π = p overlap w pred /w obs , which is always positive and less than one. 9 5 Results for aCRC data In this section we present results of applying all models (BRMA, 1a-d and 2a-d) to the illustrative example in aCRC, introduced in Section 2. Table 2 shows the between-studies correlations obtained from all the models applied to the aCRC data. The correlations are the parameters of the primary interest as they tell us about the strength of the surrogate relationships (see Sections 4.1.5 and 4.2.5). Additional results, the average effects and heterogeneity parameters, are included in Appendix D. When applying BRMA to data from all studies (regardless of the treatment contrast), the between-studies correlation is negative, -0.73 (95% CrI: -0.89, -0.49), indicating that treatment increasing the odds of tumour response leads to reduced progression rate. As shown in the scatter plot in Figure 1, the distributions of the treatment effects on the two outcomes: log OR of TR (surrogate endpoint) and log HR of PFS (final outcome), representing their association patterns, vary across the treatment contrasts. When applying NMA models 1a-c and 2a-c to the data, the between-studies correlations differ across treatment contrasts, with the highest correlation obtained for treatment contrast AC (the results are presented only for those contrasts for which at least four studies were available, the remaining results are included in Appendix D). When applying model 1b, assuming second order consistency and hence imposing additional constraints on the between-studies variance-covariance matrices, the correlations were estimated with higher precision compared to those obtained from model 1a. Applying model 1c slightly inflated the credible intervals of the correlations for some of the treatment contrasts (AC and BC) and reduced them for others (AB and BD) (compared to model 1b), which is likely due to the assumption of similar variances for the effects in each treatment arm not being reasonable for these data. The DIC value, shown in the right-hand-side column of the table, is also relatively high for model 1c and indicates poor model fit. A similar pattern is observed across models 2a-c. When assuming homogeneous variance-covariance matrices in models 1d and 2d, the between-studies correlations are high and obtained with higher precision compared to the correlation obtained from BRMA. This is likely due to the borrowing of strength across treatment contrasts, by adding information from the indirect comparisons. Based on the DIC criteria, models 1d, 2a, 2b and 2d appear to fit the data best, although models 1a and 1b appear to have a comparable fit. Prediction of the treatment effect on PFS (the final outcome) from the treatment effect on TR (the surrogate endpoint) for the treatment contrast AE, for which only one study was available (new treatment scenario), was not possible for models 1a-d, due to lack of data to estimate the surrogate relationship. Models 2a-d, assuming exchangeability of the average effects in each treatment arm across treatment contrasts, allowed us to make such a prediction. A forest plot of the observed and predicted (from BRMA and model 2d) effects on PFS is included in Figure 3. The advantage of using bvNMA methods varied across the treatment contrasts, depending on the surrogacy relationships for each contrast. Table 2 also shows the across-treatment correlation ρ t indicating, consistently across models 2a-d, a weak surrogate relationship across the treatment contrasts. Table 3 shows statistics for model comparison, introduced in Section 4.3.1, for each treatment contrast where at least four studies were available and across all the studies. For example, for the treatment contrast AB, the use of bvNMA improved the precision of the predictions in terms of the point estimate, reducing the bias from 0.23 (obtained from BRMA) to between 0.16 and 0.18. All bvNMA methods gave reduced predicted intervals for the treatment effect on PFS (from the effect on TR) compared to those obtained from BRMA, by between 11.2% and 16.5%. For contrasts AC and BC, model 1a did not contribute to reduced uncertainty of predictions, but additional assumptions of second order consistency (models 1b and 2b) and additional borrowing of strength (models 1c and 2c) or assumption of homogeneity of the correlations and heterogeneity parameters (models 1d and 2d) led to improved precision. Results from a sensitivity analysis of the data with outlying observations removed are presented in Appendix E. For the treatment contrast BD, both sets of models 1a-b and 2a-b resulted in inflated predicted intervals. The final part of Table 3 lists the average statistics across all treatments.  Table 2: Between-studies correlations corresponding to the within-treatment surrogate relationship for each model, ρ t -across-treatment correlations obtained from the models allowing for exchangeability, and DIC values corresponding to each model fitted to aCRC data. Where only one value is given for the between-studies correlation within a treatment contrast (models BRMA, 1d and 2d), the parameters are common across the treatment contrasts. A -chemotherapy alone, B -anti-VEGF therapies + chemotherapy, C -EGFRi + chemotherapy, D -EGFRi + anti-VEGF therapies + chemotherapy   Table 3: Comparison of models based on aCRC data by treatment contrast. 13 6 Illustration using simulated data 6

.1 Data simulation
To demonstrate scenarios where use of bvNMA methods has an advantage over the standard surrogacy models, data were simulated under different assumptions. In particular, we simulated data where the surrogate pattern across all studies and treatments differed from the patterns within treatment contrasts, which is detectable by mvNMA but not by BRMA. The treatment effects on two outcomes were simulated from a bivariate normal distribution: , as in model 1a. Two sets of network data were generated, each consisting of 30 studies, three treatments and three treatment contrasts with 10 studies per contrast (AB, BC and AC), under different scenarios (illustrated in Figure 4). Scenario 1 was simulated assuming weak surrogacy when ignoring treatment contrasts but strong surrogacy within each treatment contrast, with the following parameters:  Placing additional constraints on the covariance matrix by assuming second order consistency in models 1b and 1c reduced uncertainty around the correlations. Model 1d resulted in a common correlation obtained with the highest precision, however it did not take into account the differences in the between-studies variances across the treatment contrasts, in contrast to models 1a-c as shown in Appendix F. The right-hand-side column of Table 12 shows the across-treatment correlations ρ t indicating a weak across-treatment surrogate relationship. This is due to the differences in the surrogacy patterns across the treatment contrasts as well as the small number of treatment contrasts. Table 5 shows a range of statistics (described in Section 4.3.1) comparing the models in terms of their value in predicting the treatment effect on the final outcome from the treatment effect measured on the surrogate endpoint in a cross-validation procedure. The large width of the predicted interval obtained from BRMA compared to the width of the CI of the observed treatment effect estimate is due to high uncertainty, but the ratio w pred /w obs is reduced when using NMA models. Predicted intervals obtained from NMA models are between 47% and 57% narrower compared to those obtained from BRMA. The distance between the point estimate of the observed effect from the predicted effect is also much reduced when using NMA models compared to BRMA.   Table 5: Comparison of models based on simulation scenario 1.  with the variance-covariance matrix varying across treatment contrasts models the data in more detail and reveals no correlation between the treatment effects on the two outcomes within the BC treatment contrast. Figure 6 shows the predicted effects obtained from the cross-validation using (a) BRMA and (b) NMA model 1a. When using the NMA model, predictions are obtained with higher precision for contrasts AB and AC, but not BC where there was no association between the effects on the two outcomes and which is reflected by the wide predicted intervals. The across-treatment correlations ρ t in the right-hand-side column of Table 13 are obtained with high uncertainty due to the small number of treatment contrasts to estimate the correlation. Table 7 shows the statistics for the model comparison in terms of their predictive value, obtained from the cross-validation procedure. Similarly as in scenario 1, the large ratio, w pred /w obs , comparing the width of the predicted interval obtained from BRMA with the width of the CI of the observed treatment effect estimate is reduced when using the NMA models. Predicted intervals obtained from NMA models are between 19.6% and 29.3% narrower compared to those obtained from BRMA. The distance between the point estimate of the observed effect and the predicted effect is slightly reduced when using NMA models 1a-c and 2a-c compared to BRMA. These improvements, on average, are not as strong as in scenario 1, due to poor association for the treatment contrast BC. When investigating these statistics within the treatment contrasts, the improvement is largest for contrast AC where the correlation was highest (these results are presented in Appendix G).

Discussion
We have developed bivariate network meta-analytic models for surrogate endpoint evaluation to allow for a more detailed modelling of surrogate relationships within and across treatment contrasts. This methodology can help to disentangle information about surrogate relationships in data scenarios where such relationships vary across treatment contrasts and are distinct in comparison with an association pattern across treatments. Two types of surrogacy have been described by the models: surrogacy across patient populations for a given treatment contrast and surrogacy across treatments. The models will allow analysts to make predictions of the treatment effect on the final clinical outcome from the observed effect on a surrogate endpoint in a new study investigating the effectiveness of an existing treatment in a new setting or a new treatment.
There are some limitations to the models presented here. For simplicity, we focused on the models for data from two-arm studies. The methods can be extended to model multi-arm trial data in a similar manner as in Achana et al. [1]. Another limitation is related to the prior distribution for the ancillary variance-covariance matrix in model 1b, introduced to allow for the second order consistency constraints. It was based on the Cholesky decomposition which results in the prior distributions for the between-studies correlations being dependent on the ordering of the treatments in the network. We carried out a sensitivity analysis by changing the ordering of treatments in the illustrative example in aCRC; the results remained very similar to those obtained from the main analysis. Scarcity of data may also present a problem in fitting the models. For estimation of the surrogate relationships within the treatment contrasts, a relatively large number of trials per treatment contrast is needed, whereas for the across-treatment surrogacy, a range of treatments needs to be included in the network.
Nevertheless, we believe that the models have great potential for making improvements in the research area of surrogate endpoint evaluation. In our example in aCRC, the models allowed us to disentangle information on a relatively strong surrogate relationship between treatment effects on TR and PFS for the treatment contrast of EGFRi with chemotherapy vs. chemotherapy alone from a set of treatments with suboptimal overall surrogacy relationship. Moreover, in medical decisionmaking, where multiple comparisons of new health technologies against different comparators play an important role, NMA is a valuable tool in obtaining average effects across all treatment contrasts in the network of treatments. In a similar way, our proposed methodology can be used to predict the effect of a new treatment on the final clinical outcome against any comparator in the network. In conclusion, we developed a new meta-analytic method for surrogate endpoint evaluation that allows modelling of surrogate relationships in greater detail.
A Beta distribution based prior for the correlation A beta distribution was used to construct prior distributions for the between-studies correlations (all models) and for the correlation between effects on treatment arms (models 2a-c). A random variable drawn from a beta distribution r ∼ Beta(1.5, 1.5) is limited to values between 0 and 1 with probability density zero on the edges and mean value of 0.5, as seen in Figure 7 (left). Transforming this variable, such as ρ = r * 2 − 1, gives a distribution bounded by −1 and 1 with mean at zero, as shown in Figure 7 (right). This can be used as a prior distribution for the between-studies correlation, as in Burke et al (2016). The resulting prior distribution for the correlation, such as

B Surrogacy criteria
Daniels and Hughes defined the surrogacy criteria for a Bayesian meta-analytic model where the relationship between the true treatment effects on final clinical outcome µ 2i and the effect on the surrogate endpoint µ 1i was written in the form of a linear regression: The surrogate relationship between the two treatment effects, µ 2i and µ 1i , was perfect if the intercept λ 0 was zero, as then a zero effect on a surrogate would imply a zero effect on the final outcome, the slope λ 1 should not be zero for the association to be strong, with the conditional variance ψ 2 being zero. For the complete model see Daniels and Hughes (1997). A similar relationship and surrogacy criteria were described by Bujkiewicz et al (2015) in the framework of bivariate meta-analysis and extended by Bujkiewicz et al (2016) to multivariate meta-analysis. In the two papers the relationship between the regression parameters and the elements of the betweenstudies variance-covariance matrix was defined, similarly as in Bujkiewicz et al (2013). The derived relationships in the bivariate case are and If the surrogacy relationship is perfect, the conditional variance is zero: ψ 2 = 0 (Daniels and Hughes (1997), Bujkiewicz et al (2015)). Hence, from (23), τ 2 2 = λ 2 1 τ 2 1 which gives λ 1 = ± τ2 τ1 , and from (22) it implies that the correlation ρ = ±1. Also ρ 2 = 1, which some authors refer to as the study level adjusted R-squared

C Relationships between the models
The models reduce to the standard meta-analysis model for surrogate endpoints, such as the BRMA model, in a special case of data structure. When there are only two treatments in the network, as depicted in Figure 8, it can be shown that model 1a reduces to BRMA. Equations (3)-(4) become The index (12) denoting the two treatments does not vary across studies or contrasts and hence can be dropped, resulting in equations for BRMA -equations (3.1)-(3.2) in the main manuscript, with d j = β j and j = 1, 2. Tables 8 and 9 give additional resuts to the illustrative example in aCRC. Figure 3 shows the predicted effects obtained from BRMA and model 2d along with the observed estimates of the effects on PFS. The improvement in predictions was not substantial due to the weak association patterns between the treatment effects on the two outcomes.      Table 9: Mean effects and the between-studies correlations and variances for each model in the aCRC example (contrasts AE, AD and CF). A -chemotherapy alone, C -EGFRi + chemotherapy, D -EGFRi + anti-VEGF therapies + chemotherapy, E -anti-IGF1R , F -anti-IgG2 + chemotherapy 24 E Sensitivity analysis (aCRC example)

D Additional results: aCRC example
Sensitivity analysis was carried out investigating the effect of potentially influential observations (three studies with largest treatment effect on TR, due to no events in the control arm, were removed). Figure E shows the scatter plot. Tables 10 and 11 show the between studies correlations of the heterogeneity parameters.