Bayesian meta-analytical methods to incorporate multiple surrogate endpoints in drug development process.

A number of meta-analytical methods have been proposed that aim to evaluate surrogate endpoints. Bivariate meta-analytical methods can be used to predict the treatment effect for the final outcome from the treatment effect estimate measured on the surrogate endpoint while taking into account the uncertainty around the effect estimate for the surrogate endpoint. In this paper, extensions to multivariate models are developed aiming to include multiple surrogate endpoints with the potential benefit of reducing the uncertainty when making predictions. In this Bayesian multivariate meta-analytic framework, the between-study variability is modelled in a formulation of a product of normal univariate distributions. This formulation is particularly convenient for including multiple surrogate endpoints and flexible for modelling the outcomes which can be surrogate endpoints to the final outcome and potentially to one another. Two models are proposed, first, using an unstructured between-study covariance matrix by assuming the treatment effects on all outcomes are correlated and second, using a structured between-study covariance matrix by assuming treatment effects on some of the outcomes are conditionally independent. While the two models are developed for the summary data on a study level, the individual-level association is taken into account by the use of the Prentice's criteria (obtained from individual patient data) to inform the within study correlations in the models. The modelling techniques are investigated using an example in relapsing remitting multiple sclerosis where the disability worsening is the final outcome, while relapse rate and MRI lesions are potential surrogates to the disability progression. © 2015 The Authors. Statistics in Medicine Published by John Wiley & Sons Ltd.

Introduction

Surrogate endpoints are increasingly being investigated as candidate endpoints in randomised controlled trials where measuring a primary outcome of interest may be too costly, too difficult to measure or require long follow‐up time. Prior to the use of surrogate endpoint for trial design or decision‐making for regulatory or reimbursement purposes, such endpoints need to be validated. The validation takes place on three levels: by establishing a biological plausibility of the association between outcomes, assessing association between outcomes at the individual level and validating the surrogate endpoints at the study level to assess them as predictors of clinical benefit measured by the final outcome 1. It has been established that methods based on a single clinical trial are not sufficient and surrogate endpoints have to be validated based on a number of clinical trials in a meta‐analytic framework 1, 2.

A number of meta‐analytical methods have been proposed that aim to evaluate treatment effects on surrogate endpoints as predictors of the effect on a target outcome. Methods by Buyse et al. 3 were developed to model surrogate endpoints at the arm level by extending the ideas developed by Prentice 4 to a meta‐analytic framework, while those by Daniels and Hughes 2 are focussed on modelling the relationship between relative treatment effects on outcomes. The former are developed in a frequentist approach while the latter in a Bayesian framework. Various extensions to meta‐analytic approaches to evaluating surrogate endpoints have been developed, for example, for the time‐to‐event data by Burzykowki et al. 5 extended to a Bayesian framework by Renfro et al. 6.

Most methods developed to date are designed to evaluate single surrogate endpoints. In the summary of a National Institutes of Health Workshop on the use of surrogate endpoints, Gruttola et al. 7 made a number of recommendations for future research that included, for example, development of models that can accommodate measurement error, missing data and multiple surrogate endpoints and/or multiple clinical outcomes. Methods for evaluating multiple surrogate endpoints were proposed by Xu and Zeger 8 for time‐to‐event data modelled jointly with multiple biomarkers measured longitudinally but were mainly limited to individual‐level data. Other examples of validating multiple surrogate endpoints include the plasma HIV‐1 RNA and CD4+ lymphocytes as predictors of progression to AIDS in HIV‐positive patients 9, 10 and relapse rate and number of active lesions in the brain as predictors of disability progression in relapsing remitting multiple sclerosis (RRMS) 11, which used methods that did not allow for inclusion of the measurement error for the surrogate endpoint.

Bivariate meta‐analytical methods can be used to predict the treatment effect on the target outcome from the effect on surrogate endpoint (while taking into account the uncertainty around the treatment effect on surrogate endpoint ignoring of which can impact on predictions 12) as well as to combine evidence on treatment effect on both outcomes by ‘borrowing of strength’ across outcomes when evaluating new health technologies 13, 14. Extending such methods to multivariate models can be used to evaluate multiple surrogate endpoints as joint mediators of clinical benefit with a potential advantage of increasing the precision of predictions. Multivariate meta‐analysis models require within‐study correlations between multiple effect estimates which are usually not available but need to be accounted for 13, 15. If individual patient data (IPD) are not available, but results of the Prentice's criteria for surrogacy are, then we propose to use those criteria to obtain the within‐study correlations.

Models considered in this paper build on a multivariate meta‐analysis model of mixed outcomes developed by Bujkiewicz et al. 16, where the between‐study covariance is parameterised in a formulation of a product of univariate normal distributions. In this model, an assumption of conditional independence between outcomes was used to simplify the model by putting a structure on the between‐study covariance matrix. Building on this model, we extend it to two alternative models. In the first case, the assumption is relaxed to allow for full unstructured covariance to be modelled in a product normal formulation, where multiple surrogate endpoints and the final outcome are correlated and hence one surrogate endpoint can also act as a surrogate to the other. In the second case, we propose an alternative parameterisation to the one by Bujkiewicz et al. 16, which also assumes conditional independence between some of the outcomes (by putting a structure on the between‐study covariance) but is more suitable for defining criteria for surrogacy. In both models, the product normal parameterisation offers the possibility of describing the criteria for surrogacy in a greater detail compared with providing a between‐study correlation only, which would be the case when modelling the covariance structure directly. The modelling techniques are investigated using example in RRMS, where the disability worsening is the final outcome, while relapse rate and number of active MRI lesions have been considered potentially good surrogates to the disability progression 11, 17.

In the remainder of this paper, the two models (with unstructured and structured covariance matrix) are introduced in Section 2, where in addition to this a model for obtaining the within‐study correlation from published Prentice's criteria is described. In Section 2.7, the use of software is briefly reported, which is then followed by the application of the methods to data specific to the example in RRMS which is presented together with the results in Section 3. Simulation study investigating the robustness of the methods to the normality assumption as well as the performance of the methods is described in Section 4. Extensions of the two approaches to the multivariate case are developed in Section 5. The paper concludes with a discussion of methods in Section 6.

Trivariate random effects meta‐analysis with application to surrogate endpoints

Suppose in each study i, we have three estimates of treatment effect observed on each of the three outcomes Y 1, Y 2 and Y 3, where Y 3 is the treatment effect estimate for the final clinical outcome, while Y 1 and Y 2 are intermediate surrogate endpoints. Assuming the treatment effect estimates for the three outcomes in each study have a trivariate normal sampling distribution, the model can be written as In the aforementioned model, Y 1, Y 2 and Y 3 are assumed to be estimates of the correlated true treatment effects μ 1, μ 2 and μ 3 with corresponding within‐study covariance matrices Σ of the estimates (comprising of the within‐study correlations between the estimates Y and Y and within‐study standard deviations σ corresponding to each estimate Y , for each outcome j = 1,2,3 and study i). These true study‐level effects are assumed to follow a trivariate normal distribution with means (β 1,β 2,β 3) and covariance T (comprising of the between‐study correlations between the true treatment effects μ and μ and between‐study standard deviations τ corresponding to each true effect μ , j = 1,2,3 in each study i). In this hierarchical framework, Equations (1) and (2) describe the within‐study and the between‐study models, respectively.

In the remainder of this Section, the procedure for the use of the aforementioned model for the purpose of the validation of the surrogate endpoints is described in Section 2.1. The estimation of the variances for missing outcomes (omitted for the purpose of the validation) and the within‐study correlations (which are not reported for any of the studies but can be obtained from reported individual‐level surrogacy criteria) is described in Section 2.2. This is followed by setting out scenarios for modelling surrogate endpoints, by the use of an unstructured or structured between‐study covariance matrix, in Section 2.3 followed by the details of the modelling in the product normal formulation with the unstructured covariance matrix T in Section 2.4 and structured covariance in Section 2.5. The section is concluded by defining surrogacy criteria for the models in Section 2.6 and remark on the software use in Section 2.7.

Application of multivariate meta‐analysis to evaluating surrogate endpoints

Using multivariate meta‐analysis to make predictions from surrogate outcomes

Traditionally, multivariate meta‐analysis is used to estimate pooled effects for multiple outcomes from multiple studies by taking into account the correlation between the outcomes. Studies included in a multivariate meta‐analysis may all report all of the outcomes, or alternatively, it is permitted that some of the studies report only a subset of outcomes under a missing at random assumption. In the latter case, the unreported outcomes can be predicted for each of the studies by taking into account the correlation between the outcomes. In a Bayesian framework, this can be achieved by coding the unreported outcomes as missing, which are then predicted by the Markov chain Monte Carlo (MCMC) simulation of the model (1)–(2) 16. This framework can be adopted to validating surrogate endpoints (at a study level).

Using cross‐validation to examine the accuracy of predictions from multivariate meta‐analysis

The study‐level validation is conducted in the form of cross‐validation similar to the ‘leave‐one‐out’ approach as for example in Daniels and Hughes 2, except instead of taking one study out, only the estimate for the final outcome in one study is taken out. For all of the studies in the data set, one study at a time, the following procedure is conducted. For a given study, the estimate of the treatment effect on the final outcome is omitted (and coded as missing as if this particular study was a new study with the final outcome yet unknown). Then the multivariate meta‐analytic framework (1)–(2) is applied to the data including the study with the estimates of the effects on candidate surrogate endpoints known but on the final outcome missing. This missing effect is then predicted by the MCMC simulation from observed effect(s) on surrogate endpoint(s) by taking into account the data on all outcomes from the remaining studies and the relationship between the effects on all outcomes defined by the multivariate meta‐analytic model. The observed estimate is then compared with the predicted value by checking whether the value of the observed estimate falls within the predicted interval. To see if this framework will give reliable predictions for the omitted final outcome in one study, this procedure is repeated for each study i( ), giving I comparisons of how the predictions from the surrogate endpoints using the multivariate modelling approach perform. Ideally, we want to compare the performance of the model in predicting the true effect μ 3 in a new study n. But we do not know what the true effect is, so instead, we predict and then compare the estimate Y 3 with the predicted interval with the variance , where Y 1(2,3)(− denote the data from the remaining studies without the validation study n.

Post the validation process, when endpoints are established as surrogate endpoints to the final outcome, this multivariate meta‐analytic framework can be used to predict treatment effect on the final outcome from those on the surrogate endpoints in a new study, for example, for the purpose of early decision making process, in particular when the treatment effect(s) on surrogate endpoint(s) is(are) measured early compared with the final outcome.

Obtaining within‐study variances and covariances

The within‐study variability in each study i is represented by the within‐study covariance matrices Σ in model (1) which are assumed to be known. In practice, only the variances can be obtained from the reported data by calculating standard errors squared, while the within‐study correlations between the treatment effect estimates are usually not reported. Also the within‐study variance of the treatment effect for the final outcome in the study in which this effect is omitted in the cross‐validation (Section 2.1) will be unknown at the validation stage. In this case, the variance will be coded as missing and a distribution has to be placed over a missing node (required if using WinBUGS for MCMC simulation). Uniform prior distribution was placed on the missing variance, , which was recommended as a non‐informative prior distribution for variances by Lambert et al 18 (who also list other suitable prior distributions).

Within‐study covariances

Calculation of the within‐study covariances between estimates Y and Y , Σ [j,k], for each pair of outcomes j,k in each study i requires knowledge of the estimate of the within‐study correlations . These correlations between estimates of the treatment effects can be obtained by bootstrapping 2 (or double bootstrapping for correlation with uncertainty 16) of the IPD from all of the studies or a subset of studies in the meta‐analysis. Here, we consider an alternative approach where IPD is not available for any of the studies, but surrogacy on individual level has been investigated and reported by the use of Prentice's criteria.

To obtain the within‐study correlations between each pair of the treatment effect estimates Y 1, Y 2 and Y 3, calculated as a difference between two measurements in the experimental and control arms, we adopt an approach similar to the one developed by Wei and Higgins 19 who used a bivariate delta method to express the within‐study covariance between treatment effects, such as, for example, log odds ratios, in terms of the covariances between outcomes, such as the probabilities (or risks) and using a correlation between the outcomes from the literature. Here, this approach is adopted in a simpler form by representing the within‐study covariance between treatment effects in terms of the variances and correlations between the effects on specific treatment arms. In the case considered here, this is sufficient (and does not require the use of the bivariate delta method) as the correlations between effects on arm level can be obtained from Prentice's criteria of association of outcomes modelled on the absolute normal scale. However, to estimate the variances of effects in each arm, sufficient data need to be available, such as two by two table (or log odds of event with corresponding standard error) for Binomial outcomes or event rates and number of individuals (or log rates with corresponding standard errors) for Poisson outcomes.

The within‐study covariances between treatment effects on each pair of normally distributed outcomes j,k (j ≠ k;j,k = 1,2,3) can be expressed in terms of the covariances between arm‐specific effects: where Y (Y ) is the treatment difference on outcome j (k) (such as a mean difference or log odds ratio), X (X ) are the absolute treatment effects on outcome j (k) (such as mean or log odds) in arm a (where e and c stand for an experimental and control arm, respectively) in study i and where is the correlation between the absolute effects measured on outcomes j and k in each study i, which can be obtained from reported surrogacy criteria as discussed in the succeeding section.

Within‐study correlations

Taking into account the above derivation of the within‐study covariance, the within‐study correlation can be obtain from where, as noted in the previous section, it is assumed that the variances of X can be obtained from the data (for example from two by two tables for binary outcomes). Therefore, only the correlations between the effects on arm level need to be estimated in order to calculate the within‐study correlations, , between treatment effects on each pair of outcomes j and k. Consider that X is an average measurement on surrogate endpoint S , that is, X =E(S ) and X is an average measurement on final endpoint F , that is, X =E(F ), and is the correlation between those average effects. Noting that for normally distributed outcomes, the correlation between the normally distributed individual observations equals the correlation between the means ensures that should equal the correlation between the individual‐level responses for S and F that can be obtained by modelling individual patient data.

When investigating the association on individual level between surrogate endpoint S and final outcome F for patient m under treatment Z in a study i, Prentice's criteria can be used to evaluate surrogacy. The four Prentice's criteria require that (i) there is significant treatment effect on the surrogate endpoint, (ii) there is significant treatment effect on the final outcome, (iii) surrogate endpoint has a significant impact on the final outcome and (iv) the effect of treatment on the final outcome is fully mediated by the surrogate endpoint 1, 4. Similarly as by Burzykowski, Molenberghs and Buyse 1, the first two criteria for the normally distributed outcomes can be written for each study i in the following model: with correlated error structure and hence, the correlation between S and F , referred to by Buyse and Mollenberghs 20 as adjusted association (between S and F after adjustment for the treatment Z ), equals . This is the correlation in the individual‐level response for S and F which as noted above equals the correlation between the means, in this case between E(S ) and E(F ) and hence between X and X . Thus the correlation ρ equals the correlation between the absolute effects in (4) if, for example, j = S and k = F. Also the average treatment effects α and β on the two outcomes are the treatment effect estimates Y and Y in study i. The correlation ρ is required to estimate the within‐study correlation , which can be obtained from the fourth Prentice's criteria 4 which Burzykowski, Molenberghs and Buyse 1 propose to verify through the conditional distribution of the final outcome conditional on the treatment and the surrogate endpoint, where β =β −α ω /ω and γ =ω /ω . This criterion requires that if all treatment effect is mediated by the surrogate endpoint S then β should be zero. Rearranging the terms for β gives ω =ω (β −β )/α , which substituted into the formula for the adjusted association ρ gives The parameters α , β , ω and ω can be obtained by fitting model (6) to the patient data on responses S and F.

To complete the individual‐level validation, the third Prentice's criteria can be verified by with further details discussed by Burzykowski, Molenberghs and Buyse 1.

In the trivariate meta‐analysis, the within‐study correlations between treatment effects need to be estimated for each pair of outcomes. To do this using the above criteria, one of the outcomes in each pair is considered a surrogate endpoint to other outcome. Ideally, the within‐study correlations (or the Prentice's criteria from which the correlations derive) should be obtained from the IPD from each study in the meta‐analysis. When IPD is available from each of the studies, the within‐study correlation can be calculated directly. However, in the absence of the data, it may be possible to obtain the Prentice's criteria corresponding to the data from each (or some) of the studies (or some of the studies in the meta‐analysis). The Prentice's criteria provide both the assessment of surrogacy on the individual level as well as sufficient information to obtain the within‐study correlations. This approach is illustrated in Section 3, where the Prentice's criteria published in literature for one of the studies included in the meta‐analysis for the motivating example in RRMS are used to estimate the within‐study correlation. Note that establishing the within‐study correlation does not guarantee causal relationship between the treatment effects on surrogate and final outcomes. Research developments on causality of surrogate endpoints are highlighted in the Discussion section.

Scenarios for modelling of surrogate endpoints

Two scenarios for modelling the surrogate endpoints are considered here. In the first scenario, true treatment effects, μ , on all the outcomes are assumed correlated as shown in Figure 1a. In this case, the between‐study covariance matrix T in the model (2) is unstructured and can be modelled directly, by either placing inverse Wishart prior distribution on T, by Cholesky or spherical decomposition 21 (and placing separate prior distributions on the between‐study correlations and standard deviations) or by re‐parameterising the between‐study model in the product normal formulation of the series of univariate conditional distributions. The latter approach has a number of advantages. In contrast to modelling it directly by the use of Wishart prior distribution, the product normal formulation allows direct control over the prior distributions on all elements of the between‐study covariance matrix (between‐study standard deviations and correlations). It also allows to describe some association criteria which helps to assess the surrogacy in more detail compared with obtaining only between‐study correlation.

Figure 1

Scenarios for modelling surrogates endpoints: (a) all outcomes correlated giving unstructured covariance matrix T, (b) final outcome conditionally independent from the first surrogate endpoint conditional on the second giving structured covariance matrix equivalent to the precision matrix T −1 with element [1,3] equal to zero.

In the second scenario, shown in Figure 1b, in which the treatment effects on multiple outcomes may (but do not have to) be assumed to be measured sequentially in time, assumption is made that the treatment effect on the final outcome is conditionally independent from the effect on the first surrogate endpoint conditional on the effect on the second surrogate. This assumption puts a structure on the between‐study covariance matrix T resulting in element [1,3] of the precision matrix T −1 being equal to zero. This approach leads to a reduced number of parameters to estimate and is easier to implement, in particular when dealing with multiple outcomes beyond the trivariate case (for details see Sections 5.1 and 5.2).

Product normal formulation with unstructured covariance matrix

For the first scenario of the true treatment effects μ on all outcomes being correlated, represented graphically in Figure 1a, the between‐study covariance has an unstructured form. The between study model (2) for this scenario is re‐parameterised in the product normal formulation of the series of univariate conditional distributions:

Instead of placing independent non‐informative prior distributions on all the parameters and hyperparameters of the model, relationships between these parameters and the elements of the between‐study covariance matrix are derived to allow to take into account the inter‐relationship between the parameters. These relationships have the following forms: which are obtained by following the procedure described in detail in Section 5.1 for N‐dimensional case. Having established them allows to place prior distributions on the between‐study standard deviations and correlations, which is an easier task as the plausible range of values for these parameters are known or can be obtained from external sources of information (see, for example, Higgins and Whitehead 22 or Bujkiewicz et al. 16). By placing prior distributions on these parameters, for example, , τ 1(2,3)∼N(0,10)I(0,) (normal distribution truncated at value zero), the above derived relationships give the implied prior distribution on the parameters λ 21, λ 31 and λ 32 and hyper‐parameters ψ 1, ψ 2 and ψ 3. The remaining parameters are given ‘vague’ prior distributions, η 1∼N(0,1000), λ 20(30)∼N(0,1000).Note that the pooled effects β 1(2,3) on the three outcomes in model (2) are also directly linked to the model (10); β 1=η 1, β 2=λ 20+λ 21 β 1 and β 3=λ 30+λ 31 β 1+λ 32 β 2. It is possible to center the true effects on surrogate endpoints (by replacing μ with in the third and fifth line of formula (10), which can be useful if there are problems with autocorrelation) in which case the intercepts would equal the pooled effects; β 2=λ 20 and β 3=λ 30.

Product normal formulation with structured covariance matrix

A simplified model can be used by assuming conditional independence between treatment effects on some of the outcomes. For example, we can consider a situation where treatment effect is measured on different outcomes sequentially in time (or on the same outcome repeatedly in time). For example, if treatment effect on first outcome is measured at 12 months, second at 24 months and third at 36 months, then it may be conceivable to assume that the effect on the final outcome may be conditionally independent from the effect on the first outcome conditional on the second. This scenario, described in Figure 1b, leads to a simplified between‐study model with the true effect on the final outcome now conditional on the effect on the second surrogate endpoint only:

The assumption of conditional independence of the true treatment effects on outcomes three and one, μ 3 and μ 1, conditional on the effect on outcome two, μ 2 puts a structure on the covariance matrix giving the corresponding element ({1,3}) of the inverse covariance (precision) matrix equal to zero. This means that partial correlation . Because the partial correlation between μ 1 and μ 3 (adjusted for μ 2) equals , it implies that . This reduces the number of parameters in the model that need to be estimated and also simplifies the relationships between the parameters and hyperparameters of the model with the elements of the between‐study covariance matrix: which both can have an advantage in scenarios with multiple outcomes beyond trivariate case as discussed in Section 5. Similarly as in the case of the full unstructured covariance matrix, placing prior distributions directly on the between‐study standard deviations τ 1(2,3)∼N(0,10)I(0,) and correlations gives implied prior distributions placed on the parameters of the model (12), ψ 1(2,3) and λ 21(32) obtained from the the derived relationship between the two sets of parameters. The remaining parameters are given non‐informative prior distributions η 1∼N(0,1000), λ 20(30)∼N(0,1000). The pooled effects β 1(2,3) on the three outcomes in model (2) are also directly linked to the model (12); β 1=η 1, β 2=λ 20+λ 21 β 1 and β 3=λ 30+λ 32 β 2. As in the model in Section 2.4, centering of the true effects on surrogate endpoints can be applied to this models in which case the intercepts would equal the pooled effects.

Assuming such sequential structure of the treatment effects, that could be measured on the same outcome at multiple time points, the model can lead to removing some of the measurement error, for example at time two using time one. We cannot measure the true effect μ 2 at time two (only have an estimate Y 2), so having measurement at time one can improve estimate of measurement at time two. This can be useful in particular when the treatment effect at time one is measured precisely and at time two inaccurately, then the prediction at time two may be more precise (due to accounting for the correlation with outcome one) leading to better prediction at time three.

Criteria for surrogate markers

Consider first a bivariate case (with one surrogate endpoint, where the first endpoint is surrogate to the second and the third outcome is removed) as described in the first three lines of Equation (10) or (12). We can then follow the criteria set out by Daniels and Hughes 2, by which λ 21 indicates the association between the treatment effect measured by the surrogate endpoint and the treatment effect measured by the second clinical outcome (final outcome in the bivariate case with a single surrogate endpoint), therefore, we require λ 21≠0. For the association to be perfect, the conditional variance should be zero; . Also, we would expect λ 20=0 (no treatment effect on the surrogate endpoint gives no treatment effect on the target outcome), otherwise not all of the treatment effect on the target outcome is mediated by the effect on the surrogate endpoint. In the trivariate case with two surrogate markers, the effects of treatment on all biomarkers may jointly mediate the treatment effect on the final outcome. For the combined effect on the biomarkers to fully mediate the effect on the target outcome, we expect the intercept λ 30=0 and the conditional variance . The association between the effect on target outcome and each of the surrogate endpoints is expected not to be zero; λ 31,32≠0. In the sequential scenario of conditionally independent effects (described in Section 2.5), the same criteria as in the bivariate case apply, where the effect on the first outcome is a surrogate to the effect on the second and the effect on the second outcome is a surrogate to the effect on the third (final) outcome, but with additional ‘borrowing of strength’ across outcomes by taking into account the correlation structure between all of the outcomes (i.e. only conditional independence assumed).

Implementation in WinBUGS and R

All models were implemented in WinBUGS 23 where the estimates were obtained using MCMC simulation using 50000 iterations (including 20000 burn‐in). Convergence was checked by visually assessing the history, chains and autocorrelation using graphical tools in WinBUGS. All posterior estimates are presented as means with the 95% credible intervals (CrI). R was used for data manipulation and to execute WinBUGS code multiple times (for validation of surrogates for each study) using the R2WinBUGS package 24. OpenBUGS and R2OpenBUGS version of the software was used for the simulation study which was conducted using Linux (Red Hat, Inc., Raleigh, North Carolina)‐based high performance computer.

WinBUGS programs corresponding to the two TRMA models (applied to data in Tables 1 and 2) are included in Web Supplements A1 and A2.

Table 1

Data on disability progression, relapse rate and number of active MRI lesions included in the meta‐analysis.

		Disability progression					Relapse rate					MRI
Study	Follow‐up (months)	N d E	R d E	N d C	R d C		N r E	A R r E	N r C	A R r C		N m E	R m E (SE)	N m C	R m C (SE)
Paty (A)	24	124	35	124	35		124	145	124	157		124	1.80 (0.40)	124	4.9 (1.30)
Paty (B)	24	124	25	124	35		124	104	124	157		124	2.00 (0.70)	124	4.9 (1.30)
Jacobs	24	85	18	87	29		85	52	87	78		85	3.20 (0.41)	87	4.8 (0.49)
Millefiorini	24	27	2	24	9		27	12	24	31		23	3.50 (0.71)	19	7.3 (1.84)
Li (C)	24	189	57	187	69		189	172	187	240		189	9.00 (4.00)	187	15.5 (2.90)
Li (D)	24	184	50	187	69		184	159	187	240		184	5.50 (0.50)	187	15.5 (2.90)
Polman	24	627	107	315	91		627	144	315	230		627	1.90 (0.37)	315	11.0 (0.88)
Comi (E)	24	433	62	437	90		433	61	437	144		433	0.38 (0.07)	437	1.43 (0.06)
Comi (F)	24	456	69	437	90		456	68	437	144		456	0.33 (0.06)	437	1.43 (0.06)
Rudick	24	589	135	582	169		589	200	582	437		589	0.90 (0.09)	582	5.4 (0.36)
Sorensen	24	66	11	64	16		66	15	64	38		66	2.70 (0.46)	64	3.5 (0.51)
Clanet	36	400	148	402	149		400	324	402	310		400	8.00 (0.88)	402	9.0 (0.74)
Mikol	24	386	45	378	33		386	116	378	110		172	0.58 (0.11)	178	0.77 (0.18)

N d (N d ), total number of patients in experimental (control) arm with disability status recorded.

R d (R d ), number of patients in experimental (control) arm who progressed.

N r (N r ), total number of patients in experimental (control) arm with number of relapses recorded.

A R r (A R r ), number of relapses per year in experimental (control) arm.

N m (N m ), total number of patients in experimental (control) arm with number of active MRI lesions in experimental (control) arm.

R m (R m ), mean number (standard error) of active MRI lesions in experimental (control) arm.

Table 2

Indicators of individual‐level surrogacy for number of active MRI lesions and relapse rate as surrogate endpoints to disablity progression (reproduced from Sormani et al 11) and number of active MRI lesions as a surrogate endpoint to relapse rate (reproduced from Sormani et al 27).

		Prentice's criteria
Surrogate endpoint	Final outcome	1st criteriona	2nd criterionb	4th criterionc
active T2 lesions	disability progression	α 1=−0.93(0.12)	β 1=−0.37(0.19)	β S1=−0.14(0.19)
relapses	disability progression	α 2=−0.44(0.08)	β 2=−0.37(0.19)	β S2=−0.15(0.20)
active T2 lesions	relapses	α 3=−0.90(0.13)	β 3=−0.36(0.08)	β S3=−0.17(0.09)

Coefficients are reported with standard errors.

∗1st Prentice's criterion, treatment is effective on surrogate endpoint.

2nd Prentice's criterion, treatment is effective on final clinical outcome.

4th Prentice's criterion, treatment effect on final clinical outcome fully mediated by surrogate.

α 1, treatment effect on log number of MRI lesions over 1 year; α 2, treatment effect on log relapse rate

over 1 year; α 3, treatment effect on log number of MRI lesions over 2 years; β 1,β 2, treatment effect on

log odds disability progression over 2 years; β 3, treatment effect on log relapse rate over 2 years.

β and β , treatment effect on log odds disability progression over 2 years adjusted for treatment effect on

log number of MRI lesions and log relapse rate, respectively; β , treatment effect on log relapse rate

over 2 years adjusted for treatment effect on log number of MRI lesions.

Application: relapsing remitting multiple sclerosis

To illustrate the application of the methods described in Section 2, the models developed there are applied to the motivating example in RRMS described in Section 3.1. Section 3.2 (along with the Appendix A) contains details of how data in Table 1 are used to populate the models for outcomes specific to this motivating example. Results are presented in Section 3.3.

Introduction to the motivating example

To illustrate the use of the modelling techniques, we applied them to an example in RRMS. Multiple sclerosis (MS) is an inflammatory disease of the brain and spinal cord. RRMS is a most common type of MS. During the course of the disease, patients experience a series of periods of exacerbations (relapses) and remission. A large proportion of patients (25%) eventually progresses to secondary progressive disease 25. The disability progression is considered the final outcome, whereas the number of active (new or enlarging) T2 lesions in the brain obtained from the magnetic resonance imaging (MRI) and the annualised relapse rate are two potential surrogate endpoints. This example is based on work by Sormani et al. 17 who used meta‐analytic approach to evaluate whether the effects on relapse rate and number of active MRI lesions are good predictors of the treatment effects on disability progression (one surrogate endpoint at a time). In another paper, Sormani et al. 26 investigated estimates of the treatment effect on number of active MRI lesions as predictors of the effects on relapse rate. Subsequently, Sormani and colleagues used IPD to investigate individual‐level association between outcomes, where they validated the number of MRI lesions as a surrogate to the number of relapses 27. In another paper, Sormani et al. also used IPD to validate both the number of active MRI lesions and the number of relapses as surrogate endpoints to the disability progression, as individual surrogates as well as joint mediators of the treatment effect on progression 11. In the papers mentioned previously, the meta‐analytic work on study level was conducted using weighted linear regression 17, 26, whereas the association on the individual level was conducted using Prentice's criteria 11, 27.

This example with two surrogate endpoints, where one of the candidate surrogate endpoints (to the final outcome) is also a potential surrogate endpoint to the second surrogate, serves as a desirable illustration of modelling techniques investigated in this paper. We collect data on the treatment effects on those outcomes from studies included in the analysis by Sormani et al. 17. To investigate the surrogate endpoints jointly, we chose to include only the studies reporting treatment effect estimates on all of the three outcomes and those that reported sufficient information to obtain uncertainty estimates around the observed treatment effect. As a result, we obtained data from 13 studies, 11 placebo‐controlled trials and two active‐controlled trials, which are listed in Table 1 and displayed graphically in Figure 2. To obtain the within‐study correlations between the treatment effect estimates for the outcomes, we use the estimates of Prentice's criteria, reported by Sormani et al. 27 for the association between the number of MRI lesions and the number of relapses who adopted the model (6) to those outcomes on log scale and for the association of both the number of MRI lesions and the relapses with the disability progression reported by Sormani et al., 11, also on log scale. Relevant estimates are listed in Table 2. Because of limitation of the data, in this paper, the correlation is assumed constant across studies.

Figure 2

Graphical representation of data for treatment effects on MRI (surrogate endpoint 1), relapse rate (surrogate endpoint 2) and disability progression (final clinical outcome).

Scale of the outcomes and the within‐study model

The relative treatment effects on each outcome (MRI, relapse rate and disability progression) and the within‐study variances for the treatment effects in each study are calculated using data listed (and notation described) in Table 1.

The MRI effect is modelled on the log rate ratio (RR) scale, Y 1= logM R I R R= log(R m /R m ). The relative treatment effect on relapses is modelled by the log annualised relapse rate ratio (ARRR) scale, . The relative treatment effect on disability progression is modelled on the log odds ratio (OR) scale, .

The corresponding variances are obtained from the data by the use of delta method for the effects on MRI and relapses and using the standard formulae for the standard error of logO R, as shown in Appendix A.1. The within‐study correlations are obtained following methods described in Section 2.2, with details of the algebra for the data in RRMS in Appendix A.2. The correlations were obtained for one study and assumed the same in all the studies.

Note that some authors used relative risk to quantify the treatment effect on the disability progression. However, log OR scale is applied here in order to combine the summary data on effectiveness with reported Prentice's criteria, which were modelled for the progression on the log odds scale. Typically, progression is reported cross‐sectionally as a number (or proportion) of patients who progressed (i.e. whose expanded disability status scale score had increased by at least one point at the follow up). As such, the reported outcome is not suitable to model as time to event outcome on the hazard ratio scale.

Results

Applying the aformentioned formulae to the RRMS data in Table 1 and combining it with data of individual‐level surrogacy criteria in Table 2 gives the within‐study correlations: , , . The data are applied using the above models to validate surrogate endpoints as predictors of the clinical benefit on the final outcome using a cross‐validation procedure as described in Section 2.1. First, the models (1) and (10), assuming that all true effects are correlated, is applied. Results of this application are presented in Section 3.3.1. Then the models (1) and (12), assuming the conditional independence between the effects on the final outcome and the first surrogate endpoint is applied and results presented in Section 3.3.2.

Results from the model with unstructured between‐study covariance matrix

To quantify the surrogacy criteria overall, the model was first applied to the full data (no outcome missing for any of the studies). The results (parameters with the 95% CrIs) are shown in Table 3 (top). The coefficients λ 20 and λ 21 and variance represent the association between the treatment effect on the number of active MRI lesions and the effect on the relapse rate, whereas the coefficients λ 30, λ 31, λ 32 and describe how the effects on number of MRI lesions and relapse jointly mediate the effect on the disability progression (as described in Section 2.6). The parameters on the left‐hand‐side indicate that there was a poor association between the effects on MRI and effects on relapses, as the interval of the slope contains zero and the variance appears greater than zero. The parameters on the right‐hand‐side suggest that the association between the effects on both surrogate endpoints (MRI and relapse rates) are not strongly associated with the effect on disability progression (slopes contain zero), however, the variance appears small suggesting that the effect on MRI and relapse rate largely mediate the effect of the treatment on disability progression. In both cases, the CrIs for the intercepts contain zero, which is encouraging in the sense that zero treatment effect on MRI would be expected to predict zero effect on relapse rate in the first case and the zero effects measured by the MRI and relapse rate would lead to predicted zero effect on the disability progression in the second case. However, in both cases the CrIs are wide, hence, the predictions may not be accurate.

Table 3

Surrogacy criteria obtained from the trivariate models applied to the full data (no missing outcomes).

Parameter	Mean (95% CrI)	Parameter	Mean (95% CrI)
Unstructured between‐study covariance
λ 20	−0.28 (−0.73, 0.16)	λ 30	−0.09 (−0.38, 0.21)
λ 21	0.26 (−0.12, 0.65)	λ 31	0.00 (−0.38, 0.31)
		λ 32	0.43 (−0.02, 1.06)
ψ22	0.15 (0.05, 0.37)	ψ32	0.02 (0.00, 0.10)
Structured between‐study covariance
λ 20	−0.24 (−0.70, 0.20)	λ 30	−0.06 (−0.31, 0.19)
λ 21	0.30 (−0.09, 0.70)	λ 32	0.48 (0.11, 0.88)
ψ22	0.15 (0.05, 0.36)	ψ32	0.02 (0.00, 0.10)

After estimating the surrogacy criteria using the full data, the cross‐validation process was carried out. The results of the validation are presented in Table 4 (column four, along with the results from the bivariate model and trivariate with structured covariance matrix), which shows the values of observed effects on the final outcome in each study (column two) and the effects on the final outcome predicted by the model from the effects on the surrogates (given full data on all other studies) with the 95% CrIs. In each row, the values correspond to the validation of prediction for a study whose corresponding observed effect on the final outcome is in the second column. In the take‐one‐out approach of the cross‐validation procedure, the regression parameters (the intercepts, slopes and conditional variances, as in the between‐study model (10)) did not change substantially. The full set of those values is included in the Web Supplement B. The predicted intervals contain the observed value of the treatment effect on the final outcome for all of the studies confirming good fit of the model.

Table 4

Comparison of results of validation obtained from the two trivariate models and a bivariate model.

		Predicted
Study	Observed	Bivariate	Trivariate Unstructured	% red.	Trivariate Structured	% red.
Paty (A)	1.00 (0.57, 1.74)	0.86 (0.42,1.74)	0.86 (0.43, 1.73)	0.67	0.88 (0.44, 1.73)	3.16
Paty (B)	0.64 (0.36, 1.16)	0.76 (0.37, 1.56)	0.77 (0.38, 1.54)	2.24	0.77 (0.38, 1.54)	3.29
Simon	0.53 (0.27, 1.06)	0.78 (0.36, 1.73)	0.80 (0.36, 1.77)	0.19	0.79 (0.37, 1.72)	2.10
Millefiorini	0.13 (0.02, 0.67)	0.62 (0.11, 3.46)	0.65 (0.12, 3.62)	−0.29	0.62 (0.11, 3.44)	0.32
Li (C)	0.73 (0.47, 1.12)	0.79 (0.43, 1.45)	0.79 (0.45, 1.41)	4.46	0.80 (0.45, 1.42)	5.19
Li (D)	0.63 (0.41, 0.98)	0.79 (0.43, 1.44)	0.81 (0.44, 1.50)	−1.07	0.79 (0.45, 1.39)	7.23
Polman	0.50 (0.36, 0.69)	0.59 (0.32, 1.11)	0.61 (0.35, 1.07)	9.67	0.58 (0.33, 1.03)	8.62
Comi (E)	0.64 (0.45, 0.92)	0.62 (0.35, 1.11)	0.62 (0.36, 1.08)	3.71	0.62 (0.36, 1.07)	5.13
Comi (F)	0.69 (0.49, 0.97)	0.64 (0.36, 1.13)	0.64 (0.38, 1.09)	5.77	0.63 (0.37, 1.08)	5.93
Rudick	0.73 (0.56, 0.95)	0.62 (0.37, 1.05)	0.61 (0.37, 1.01)	4.18	0.61 (0.38, 0.98)	10.01
Sorensen	0.57 (0.24, 1.36)	0.62 (0.24, 1.64)	0.64 (0.23, 1.82)	−7.52	0.63 (0.24, 1.66)	0.63
Clanet	1.00 (0.75, 1.33)	0.89 (0.49, 1.63)	0.88 (0.49, 1.58)	3.21	0.92 (0.53, 1.60)	8.30
Mikol	1.39 (0.87, 2.23)	0.82 (0.44, 1.55)	0.83 (0.45, 1.53)	2.65	0.85 (0.46, 1.57)	3.45
average % reduction in CrI				2.14%		4.87%
DIC			350.1		347.9

The % red refers to the percentage reduction in the width of the credible interval corresponding to the prediction from the trivariate model, with the unstructured (columns 4 and 5) or structured (columns 6 and 7) between‐study covariance matrix, compared with the width of the interval corresponding to the prediction from the bivariate model (column 3).

DIC, deviance information criteria.

Results from the model with structured between‐study covariance matrix

Models (1) and (12) with structured between‐study covariance matrix, where the effect on the disability progression (the final outcome) is conditionally independent from the effect on the MRI (the first surrogate endpoint) conditional on the effect on relapse rate (the second surrogate endpoint) is now used for validation. Similarly as in the case of the model with unstructured covariance matrix, the surrogacy criteria are estimated using the full data, which is then followed by the cross‐validation. In Table 3 (bottom) the parameters are shown together with the 95% CrIs for the model applied to the full data. The parameters λ 20, λ 21 and describe the association between the treatment effects on the MRI and the relapse rate, while λ 30, λ 31 and the association between the treatment effects on the relapse rate and the disability progression (conditional on the effect on MRI). The association between the effects on the MRI and the relapse rate is not strong which is indicated by the interval of λ 21 containing zero and the variance significantly larger than zero. However, the association between the effect on the relapse rate, as a surrogate endpoint, with the effect on the disability progression, as the final outcome, (conditional on the effect on MRI) appears strong as indicated by the non‐zero slope λ 31 and the small variance . Also the interval of the intercept λ 30 containing zero indicates that zero effect on the relapse rate is likely to imply zero effect on the disability progression. In the take‐one‐out approach of the cross‐valdation procedure, the regression parameters (the intercepts, slopes and conditional variances, as in the between‐study model (12)) did not change substantially. The full set of those values is included in the Web Supplement B. Table 4 shows the results of cross‐validation (with column six corresponding to predictions from the model with structured between‐study covariance). The results of predictions are similar to those obtained from the model with unstructured between‐study covariance matrix.

Comparison of the results from the two models and those from a bivariate model

When carrying out the cross‐validation process, we want to ensure that not only predicted CrIs contain the actual observed values but also that the intervals are narrow. Inclusion of multiple surrogate endpoints can potentially lead to reduced intervals and hence better predictions. To compare the aforementioned results of the validation from the two trivariate models, they are shown side by side in Table 4 as well as graphically in Figure 3 alongside those from a bivariate model with the effect on the relapse rate as a surrogate for the effect on the disability progression. In Table 4, apart from the predicted values, the percentage reduction in the width of the credible interval relative to the width of the interval obtained from the bivariate model is shown for both trivariate models. On average, the model with the unstructured between‐study covariance matrix gave intervals 2% narrower than those from the model using a single surrogate endpoint, whereas the model with the structured covariance matrix led to 5% reduction in uncertainty. Although the average gain in precision is modest, the inclusion of both surrogates does improve the predictions and this reduction in uncertainty may be sufficient to improve the decision making process based on such predictions. This may be the case in particular for larger studies, such as by Polman, Rudick and Clanet in our illustrative example, where the gain in precision was also larger, up to 10%, likely due to the treatment effect on the surrogate endpoints being measured with larger precision compared with the smaller studies. Inclusion of multiple surrogate endpoints may lead to a more substantial gain in precision in other disease areas and when data on some outcomes in some of the studies are missing at random 28. Deviance information criteria (DIC) obtained for the two trivariate models using the complete data, showed at the bottom of Table 4, suggests that the fit by both models is comparable. The model assuming conditional independence is simpler and easier to implement as shown in Section 5. It also requires fewer parameters to estimate as discussed further in Section 6. If such assumption can be justified, the model with the structured covariance may be more practical.

Figure 3

Forest plot, showing for each study the observed value of the OR of disability progression with corresponding confidence interval (CI) and the predicted values with corresponding credible intervals (CrIs) from a bivariate model, trivariate model with unstructured covariance matrix and from the trivariate model with structured covariance matrix.

Simulation and sensitivity analysis

Methods

A simulation study was carried out to compare the performance of the three models: a bivariate model with a single surrogate endpoint and the two trivariate models with two surrogate endpoints, the model with unstructured between‐study covariance matrix and with the structured covariance matrix. Data were simulated under two scenarios considering alternative covariance structures: unstructured, by simulating data from the model (1)+(10) and structured by simulating data from the model (1)+(12). The parameters of the between‐study model were set to be comparable with those corresponding to the RRMS data, namely, τ 1,2,3=0.5 and and η 1=−0.3.

The regression coefficients and the conditional variances were obtained from the between‐study standard deviations and correlations following the formulae (11) and (13) for the models with unstructured and structured covariance matrices, respectively. The within‐study correlations were set to the same values as in the example in RRMS (obtained from the Prentice's criteria). The within‐study variances were simulated by drawing the precisions (inverse variances) from the gamma distributions; σ 1(2,3)=1/P 1(2,3), P 1(2,3)∼Γ(α 1(2,3),θ 1(2,3)), where α 1(2,3) are the shape parameters and θ 1(2,3) the scale parameters, which were obtained using the method of moments: E(P 1,2,3) = α 1,2,3/ξ 1,2,3, , where ξ 1,2,3=1/θ 1,2,3 is a rate parameter. By summarising the inverse variances from the RRMS data, the following parameters were obtained: E(P 1) = 30, E(P 2) = 150, E(P 3) = 25, V(P 1) = 420, V(P 2) = 15000 and V(P 2) = 275, giving the following shape and rate parameters: α 1=2.14, α 2=1.5, α 3=2.3, ξ 1=0.07, ξ 2=0.01 and ξ 3=0.09. Because of the structure of the gamma distribution, some of the simulated precisions were very close to zero, resulting in very large variances. This led to some problems with the estimation. To overcome this issue, a constraint was placed on the simulated value of the precision by discarding the precisions resulting in variances larger than 2 (this number was taken as an arbitrary cut off, large enough to be much larger than the variances in the RRMS data and hence including all plausible variances in the population but small enough not to produce problems with the estimation). R code for the data simulation is included in the Web Supplement C.

The cross‐validation procedure is applied to each simulated data set. This time, however, the true effect μ 3 in a validation study n is known as it has been simulated, so the cross‐validation can be performed on the true effects (which in real circumstances we would like to predict) by comparing the simulated μ 3 with predicted interval of with the corresponding variance . The two trivariate models and also the bivariate model are applied to the simulated data sets to compare their performance, by estimating bias of the mean , root‐mean‐square error (RMSE), coverage of 95% credible intervals (CrI) and the potential reduction of the width of the predicted interval by calculating the ratios of the width of the interval from a trivariate model w 3 (with two surrogate endpoints) with the width of the predicted interval from the bivariate model w 2.

Sensitivity analysis

To investigate the sensitivity of the methods to the assumption of normality of the data, another simulation study was carried out where the data were simulated from multivariate t‐distribution as well as mixture normal distribution.

Simulation study with t‐distribution

To simulate data from thet‐distribution with unstructured covariance matrix, the true treatment effects on the three outcomes were generated from the following model: and from the t‐distribution with structured covariance matrix, from the following model where and d f = 4. The individual study estimates were simulated from the trivariate t‐distribution The simulation was conducted in R software using the rt command (and scaling the simulated data by ) for univariatet‐distributions in the between‐study models (14) and (15) and the rmvt in the within‐study model (16).

Simulation study with mixture normal distribution

Data with more severe departure from normality were generated by the use of mixture normal distributions. This was achieved by replacing the univariate normal distribution for the true treatment effect on the first outcome in the models (10) and (12) with the mixture normal distribution with p 1=0.5, p 2=0.3 and p 3=0.2. This deviation from normality feeds through to the true effects μ 2 and μ 3 by the linear association of those effects with μ 1. These now non‐normal true effects μ 1(2,3) are then used as mean values when generating the within study data from the multivariate normal distribution giving data with ‘distorted’ normality.

Data sets including the treatments effects (and corresponding sampling variances) on three outcomes in 15 studies were generated in 1000 simulations for each scenario. 0.1% of simulation runs were discarded because of precisions resulting in too high variances, as explained in the methods Section 4.1. Simulation results are presented in Table 5. The bias of mean predicted effect for a validation study n was comparable across all models and data scenarios. The RMSE was larger when data were simulated from a model with unstructured covariance matrix, regardless of the distribution or model fitted to the data. All models (the bivariate and the two trivariate) seemed to perform equally well, giving coverage of 95% credible interval close to 95% for most scenarios, except for the data generated from the mixture normal distribution where the coverage was slightly inflated to 98% (because of the three normal distributions being approximated by models with a single normal distribution, leading to the inflated variance of predictions). For data generated from either multivariate normal or t‐distribution with unstructured covariance matrix, both trivariate models gave on average 4% reduction of the width of the predicted interval compared with the intervals obtained from the bivariate model. When data was simulated using normal or t‐distribution with structured covariance matrix, the trivariate model with unstructured covariance gave on average 2% reduction in uncertainty around the predicted effect compared with the uncertainty around predictions obtained from the bivariate model, while the model with structured covariance matrix gave predictions with intervals 5% narrower compared with those obtained from the model with a single surrogate endpoint. When data were simulated from the mixture of normal distributions, the trivariate model with unstructured covariance did not produce any gain in precision of predictions, while the model with structured covariance matrix gave predicted intervals that were on average 9% narrower compared with those obtained from the bivariate model.

Table 5

Results of simulation studies.

Meta‐analysis model	Bias of mean μ^_3n	RMSE of μ^_3n	Coverage of 95% CrI for μ^_3n	Median w 3d/w 2d
Scenario 1: Data simulated from normal TRMA with UCM
BRMA	−0.002	0.46	0.96
TRMA UCM	0.003	0.45	0.96	0.96
TRMA SCM	−0.001	0.46	0.95	0.96
Scenario 2: Data simulated from normal TRMA with SCM
BRMA	0.006	0.36	0.95
TRMA UCM	0.003	0.36	0.95	0.98
TRMA SCM	0.006	0.35	0.94	0.95
Scenario 3: Data simulated from TRMA with UCM and t‐distribution
BRMA	−0.002	0.48	0.95
TRMA UCM	−0.004	0.47	0.95	0.96
TRMA SCM	−0.002	0.48	0.94	0.96
Scenario 4: Data simulated from TRMA with SCM and t‐distribution
BRMA	−0.0001	0.36	0.95
TRMA UCM	0.002	0.37	0.95	0.98
TRMA SCM	0.0002	0.36	0.94	0.95
Scenario 5: Data simulated from TRMA with UCM and mixture normal
BRMA	−0.007	0.49	0.98
TRMA UCM	0.003	0.47	0.98	1.00
TRMA SCM	−0.001	0.47	0.97	0.91
Scenario 6: Data simulated from TRMA with SCM and mixture normal
BRMA	−0.002	0.36	0.98
TRMA UCM	0.0001	0.37	0.98	1.01
TRMA SCM	−0.0001	0.36	0.97	0.91

RMSE, root‐mean‐squared‐error; CrI, credible interval

UCM, unstructured covariance matrix; SCM, structured covariance matrix; TRMA, trivariate random effects meta‐analysis;

w 3 (w 2), width of the predicted interval from TRMA (BRMA)

MC errors of the predicted mean effects were less than 0.012 in scenarios 1–4 and less than 0.015 and 0.025 in scenarios 5 and 6, respectively

Discussion of the results

The outcomes of the simulation study were broadly in agreement with those obtained from the case study. The predicted intervals obtained from TRMA models were narrower compared with those obtained from BRMA, but this reduction was less pronounced when using the TRMA UCM model on the data simulated from a model with structured covariance matrix (one of the surrogates correlated to the other surrogate but less so to the final outcome; scenarios 2, 4 and 6 in Table 5 and also in the RRMS case study). Using TRMA UCM on data simulated from the same model gives the same reduction in uncertainty as when using TRMA SCM but this effect of the addition of the second surrogate on the uncertainty around the predicted effects diminished with the departure from normality (scenario 5). There was no effect of adding the second surrogate endpoint on the RMSE which was almost the same across the three methods within each scenario. When data were simulated from non‐normal distribution (t‐distribution or mixture normal) with structured covariance matrix, the RMSE was slightly larger when using TRMA UCM (compared with BRMA or TRMA SCM). This is likely due to the TRMA UCM model forcing too rigid correlation structure on the data leading to bias when making predictions for outlying observations in data with not as strong correlation pattern. When data represents all outcomes correlated (from a distribution with unstructured covariance matrix), both TRMA models seem to perform equally well, with slightly smaller RMSE when using TRMA UCM if the data are normally distributed. However, for non‐normally distributed data, gain in precision is only present when using TRMA SCM. When data corresponds to the scenario with the structured covariance matrix the TRMA SCM model seems to perform better than TRMA UCM in terms of both RMSE and the uncertainty of the predictions.

Multivariate random effects meta‐analysis

Methodology introduced in Section 2 can be extended to a scenario with multiple surrogate endpoints. Suppose, we have estimates of treatment effects observed on N outcomes, in each study i, and Y is the final clinical outcome of interest, while are intermediate surrogate endpoints. If the estimates of the treatment effects on all of the outcomes are assumed normally distributed and correlated, then they follow a multivariate normal distribution: In the aforementioned model, outcomes Y 1, Y 2, and Y are assumed to be estimates of correlated true effects μ 1, μ 2, and μ with corresponding within‐study covariance matrices Σ of the estimates. These study‐level effects follow a multivariate normal distribution with means and covariance T, In this hierarchical framework, Equations (18) and (19) describe the within‐study and the between‐study models, respectively. The between‐study model can be reparameterised to extend the scenarios for modelling of surrogate endpoints described in Sections 2.3, 2.4 and 2.5 to the multivariate case, using the unstructured and structured covariance matrices.

Assuming that true treatment effects on all the outcomes are correlated, we can parameterise the between‐study model (19) in the form of product normal formulation by extending model (10) to

In this model, prior distributions need to be placed on all the parameters. Non‐informative normal distributions are placed on the mean effect η 1∼N(0,1000) and the intercepts . Similarly as in the trivariate case, we place prior distributions on the between‐study correlations and between‐study standard deviations (elements of matrix T in (19) for which we are more likely to anticipate a range of values or, in some applications can obtain an external information to construct informative prior distributions for them). The relationships between the model hyper‐parameters (conditional variances and slopes ) and the between‐study parameters (correlations and standard deviations) give implied prior distributions for those hyper‐parameters and also ensure that the between‐study covariance is positively defined, These relationships are obtained by calculating the variances and correlations in terms of the hyperparameters and rearranging the equations for the variances and solving the set of simultaneous equations for correlations. Details are given in Appendix B.1. This becomes a complex task in higher dimensions. Alternative and simpler model can be used by assuming a structure of the between‐study covariance matrix as described in Section 5.2.

Criteria for surrogate markers

In the case of multiple endpoints with N − 1 surrogate markers, the effects of treatment on all biomarkers may jointly mediate the treatment effect on the target outcome. For the combined effect on the biomarkers to fully mediate the effect on the target outcome, we expect the intercept λ =0 and the conditional variance . The association between the effect on target outcome and each of the surrogate endpoints is expected not to be zero; λ ≠0, ( ). Similarly, the relationship between effects on other outcomes can be investigated. For example, if we consider outcomes Y 1 to Y to be potential surrogates to outcome Y (M < N), then coefficients for η , namely, and , are investigated in the same manner, giving an option for considering multiple ‘final’ endpoints at the same time.

Scenario with two surrogate endpoints and final outcome measured ‘in sequence’, described in Section 2.5 can be extended to the multivariate case. If we imagine that the N outcomes are ordered in a sequence (for example according to measurement time or other reasons that would impose such correlation structure), a conditional independence between any pair of outcomes that are not ‘neighbours’ can be assumed, conditional on the outcomes placed in the sequence in between that particular pair.

This leads to a structure being placed on the between‐study covariance matrix in such a way to fully take into account of the correlations between the treatment effect on pairs of outcomes (for example those that are measured one after another in a time sequence, but assume conditional independence of other effects). The elements of the precision matrix T −1 corresponding to the effects that are conditionally independent become zero and only those on diagonal and immediate off‐diagonals are non‐zero. The between‐study model (19) is then parameterised in the product normal,

Similarly as in previous models (the trivariate and the multivariate with unstructured between‐study covariance), the parameters of the above model can be expressed in terms of the elements of the between‐study covariance matrix T (19): and obtained from the formulae for the correlations in Appendix B.2. Non‐informative prior distributions are then placed on the between‐study correlations and standard deviations, , , as well as the remaining, independent parameters, η 1∼N(0,1000), . In this form, the model is much easier to implement compared with the model with the unstructured covariance matrix in Section 5.1.

Discussion

We have developed a multivariate meta‐analytic framework to include multiple surrogate endpoints when predicting the treatment effect on the final clinical outcome in an early drug evaluation process. The validation process discussed here aims to evaluate how well the effect of a treatment on multiple surrogate endpoints can jointly predict treatment effect on the final clinical outcome. Two approaches were developed, one assuming all effects being correlated giving an unstructured between‐study covariance matrix and a model assuming conditional independence between the effects on some of the endpoints giving a structured covariance matrix. The first model makes fewer assumptions about the data but requires estimation of a large number of parameters (which may be difficult if the number of studies is relatively low) and also may be difficult to implement in higher dimensions. The second model makes an assumption of conditional independence of some effects but leads to a reduced number of parameters to estimate and is easier to implement as the relationships between the parameters of the model and the elements of the between‐study covariance matrix have a simple form. For example, in a scenario with five endpoints (four surrogate endpoints and one final outcome), the model with unstructured covariance matrix is set up to estimate five between study variances and ten between‐study correlations, while the model assuming conditional independence with structured between‐study covariance matrix estimates only four between‐study correlations. However the assumption of conditional independence of some of the outcomes may lead to underestimating some correlations which may impact on predictions. This modelling framework, however, allows for flexible modelling of the correlation structure where a choice can be made about which correlations need to be fully taken into account. A model ‘in between’ those two discussed can be implemented, such as the one showed schematically in Figure 4 which extends the scenario with five endpoints in sequence (structured covariance matrix model) by taking into account of the correlation between the final and the second to last outcome (now removing the zero from the element [3,5] of the precision matrix T −1 which was present in the sequential model with conditional independence of effects on outcomes three and five). The desirable correlation is taken into account, while an assumption of conditional independence of the remaining outcomes is still made thereby reducing the number of correlations from ten in the model with unstructured covariance matrix to only five. See all three scenarios represented graphically in the Web Supplement D. The results of the case study and the simulation study suggest that the choice of the model (in terms of the covariance structure) when made in line with the correlation structure of the data will lead to better predictions.

Figure 4

Example of a scenario of modelling multiple surrogate endpoints with a choice of a correlation structure.

Further research is required to extend the aforementioned methodology to more exact methods for Binomial or count data to relax the assumption of normality, for example, by adopting methods similar to those introduced by Stijnen et al. 29 using exact likelihood, by means of a generalised linear mixed models. However, our simulation study showed robustness of the methods to the departures from normality of the data. Another limitation is the prediction of the treatment effect ‘estimate’ rather than the true effect. Ideally, we want to know how our predicted treatment effect compares with the true treatment effect. However, although the models have the ability to predict the true effects, such effect cannot be used for validation based on the real data as the true treatment effect is unknown which limits the cross‐validation to comparing the predicted treatment effect ‘estimates’ with the observed (but pretended missing) treatment effects for the final outcome.

The approach to validation of surrogate endpoints presented here focusses mainly on validating the endpoints as good (joint) predictors of the clinical benefit on the final outcome. It is now well established that the surrogacy has to be validated on both individual and study level. Prentice's criteria can be used to validate surrogate endpoints on the individual level when IPD are available. The use of the Prentice's criteria to estimate the within‐study correlations between the treatment effects on different outcomes (by assuming one endpoint in each pair of outcomes can be considered a surrogate to the other) provides a bridge between validation on the two levels. This approach may be useful when the sharing of patient data from multiple clinical trials is problematic, but owners of the data may be willing to share Prentice's criteria.

There are some limitations of Prentice's criteria for surrogacy. As discussed by a number of authors 30, 31, 32, satisfying those criteria does not guarantee a causal relationship between the treatment effects on the surrogate and the final outcome. An example of data scenario where Prentice's fourth criterion is satisfied but it does not support a claim of causality can be found in Buyse et al, who discuss a number of approached to evaluation of surrogate endpoints 32. As pointed out by Joffe and Greene 30, the statistical model describing the Prentice's criteria does not account for the common causes of the surrogate and the final outcomes. The fourth Prentice's criteria, sometimes referred to as the surrogacy criteria, is based on the assumption that the treatment effect can be partitioned into direct and indirect effects with the indirect effect being the part mediated by the surrogate and the direct being the part of the effect not mediated by the surrogate. Ignoring the common causes by modelling the direct effect on the final outcome by conditioning on the effect on the surrogate is, however, not equivalent with experimentally controlling the surrogate, as pointed out by Pearl 33. As discussed by VanderWeele 31, this can lead to the ‘surrogate paradox’ where despite the positive association between treatment and the surrogate and the surrogate and the final outcome, the association between the treatment and the final outcome may be negative. Chen et al 34 referred to this phenomena as effect reversal.

This issue of the causal effect mostly affects surrogate evaluation when such analysis is based on a single study. The meta‐analytic approach, such as described in our paper, follows the causal association paradigm which, as discussed by Joffe and Greene, is based on establishing the association between the treatment effects on the candidate surrogate endpoint and on the final outcome (rather than modelling the effect of surrogate on the final outcome) 30. The authors point out that this approach is more useful for evaluation of surrogate endpoints as it is free from the restrictions of the causal effect paradigm as the causal association does not require experimental manipulation of the surrogate. The meta‐analytic approach to the causal association study involves associations between quantities that derive directly from randomisation and as such are average causal effects. Meta‐analytic approach, such as proposed in our study, is based on data from a number of studies or subgroups and is likely to include heterogeneous treatment contrasts which is an obvious advantage over evaluation based on a single study. Such a single trial validation cannot guarantee that an association between effects confirmed based on individual data under one treatment will hold in other interventions. Moreover, Alonso et al, who investigated the relationship between the causal inference and meta‐analytic approaches to the surrogate endpoint evaluation, have shown that causal effect of the surrogate on the final outcome validated based on a single study may not be confirmed in the meta‐analytic setting, in particular, when the between‐study heterogeneity is large and the causal effect is weak 35. The authors have also concluded that a surrogate endpoint that successfully validated in a meta‐analytic setting, on both individual and trial level, is likely to be confirmed when evaluated based on causal inference framework.

In this paper, we used the Prentice's criteria to estimate the within‐study correlation between the treatment effects, as ignoring the within‐study correlation has known consequences 15. Moreover, ideally, the correlation should be obtained for each study independently to account for potential differences in the association depending on the treatment contrast. The main goal in the development of these models was building a framework under which all available evidence can be used to predict future treatment effect on the final outcome while taking into account of uncertainty around relevant parameters. In particular, these approaches have advantage over simple regression models as they take into account measurement error of the treatment effect on surrogate endpoints ignoring which can impact on predictions 12. Further research on putting this work in the framework of causal mediation could provide solution to some of the aforementioned limitations.

supporting Info Item

Click here for additional data file.