Research Article Received 2 July 2013,

Accepted 9 March 2015

Published online in Wiley Online Library

(wileyonlinelibrary.com) DOI: 10.1002/sim.6492

Network meta-analysis of longitudinal data using fractional polynomials J. P. Jansen,a,b*† M. C. Vieirac and S. Coped Network meta-analysis of randomized controlled trials (RCTs) are often based on one treatment effect measure per study. However, many studies report data at multiple time points. Furthermore, not all studies measure the outcomes at the same time points. As an alternative to a network meta-analysis based on a synthesis of the results at one time point, a network meta-analysis method is presented that allows for the simultaneous analysis of outcomes at multiple time points. The development of outcomes over time of interventions compared in an RCT is modeled with fractional polynomials, and the differences between the parameters of these polynomials within a trial are synthesized across studies with a Bayesian network meta-analysis. The proposed models are illustrated with an analysis of RCTs evaluating interventions for osteoarthritis of the knee. Fixed and random effects second order fractional polynomials were applied to the case study. Network meta-analysis with models that represent the treatment effects in terms of several parameters using fractional polynomials can be considered a useful addition to models for network meta-analysis of repeated measures previously proposed. When RCTs report treatment effects at multiple follow-up times, these models can be used to synthesize the results even if reporting times differ across the studies. Copyright © 2015 John Wiley & Sons, Ltd. Keywords:

network meta-analysis; mixed treatment comparison; repeated measures; fractional polynomials; study level data

1. Introduction Although often placed at the top of evidence hierarchies, single randomized controlled trials (RCTs) rarely provide adequate information for addressing the evidence demands of patients, clinicians, and policy makers. Instead, each trial provides a piece of evidence that, taken together, addresses important questions of healthcare decision-makers [1]. Traditional pair-wise meta-analyses of RCTs are increasingly used to synthesize the results of different trials evaluating the same intervention(s) to obtain an overall estimate of the treatment effect for the intervention relative to the control. In the last decade, network meta-analysis has been introduced as a generalization of pair-wise metaanalysis [2–10]. When the available RCTs of interest not all compare the same interventions – but each trial compares only a subset of the interventions of interest – it is possible to represent the evidence base as a network where all trials have at least one intervention in common. Such a network of trials involving treatments compared directly and indirectly can be synthesized by means of network metaanalysis [3–10]. Many RCTs report treatment effect estimates for the outcomes of interest at multiple time points. Possibilities for the meta-analysis of trials with repeated measures include analysis of one relevant time point that is common across all or most studies; analysis at the final time point of each trial; or separate meta-analysis at each relevant time point. The obvious disadvantage of these methods is that not all time points are considered simultaneously, which does not allow for the trend of treatment effects over time to be evaluated. Furthermore, choices regarding the time points to be analyzed may affect the evidence

a Redwood Outcomes, San Francisco CA, U.S.A. b Tufts University School of Medicine, Boston MA, c Previously at Mapi, Boston MA, U.S.A. d Mapi, Toronto, Canada *Correspondence to: Jeroen P. Jansen, † E-mail: [email protected]

U.S.A.

Tufts University School of Medicine, 145 Harrison Ave Boston, MA 02111, U.S.A.

Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

base used and therefore the findings.Another problem is the issue of multiple testing with multiple metaanalyses on the same data. Finally, the correlations between the time-specific treatment effects are ignored when, repeatedly, a meta-analysis is performed for each consecutive time point [11]. Hence, a metaanalysis that synthesizes the treatment effects at multiple time points simultaneously and accounts for the correlation over time is preferable. In order to estimate the correlation of observations over time within RCTs, individual patient level data (IPD) meta-analysis can be considered the gold standard. Unfortunately, researchers rarely have access to IPD for all the relevant RCTs. When IPD are not available, the ideal aggregate (AD) or study level data include the treatment effect estimates at multiple time points, along with their variance and covariance estimates [12]. However, publications of RCTs rarely report covariance estimates of treatment effects over time, which means the within-study correlation in an AD meta-analysis of repeated measures is often ignored or assumed to take on a certain value (tested through sensitivity analyses), possibly based on external evidence [12–14]. To evaluate the overall trend of treatment effects over time based on AD meta-analysis, Peters and Mengerson have illustrated that synthesizing study-specific trends based on treatment effects at each time point within studies is preferred over pooling effects across studies at each time-point, followed by fitting a trend line [14]. With the latter approach, violation of the independence assumption (i.e. ignoring within-trial correlation) will have a greater impact in terms of overly precise trend estimates. Recently, several models for network meta-analysis of repeated measures have been presented. Wandel et al. introduced a repeated measures network meta-analysis model for continuous outcomes over time. A limitation of their model is that treatment-specific outcomes at each time point for each trial are assumed to be exchangeable around a trial and treatment-specific mean result that is constant over time [15]. While this model described this particular dataset well, in other applications, there may be more structure in the repeated measurements that should be accounted for in the model. Dakin et al. presented fixed and random effects models for network meta-analysis of continuous outcome over time assuming constant or time-point specific relative treatment effects in the presence or absence of a constant constrained baseline [16]. Because Dakin et al. did not assume a functional form for the development of relative treatment effects over time (other than assuming constant treatment effects from one time point or fixed period to the next), their models cannot compare trends of treatment effects between treatments. Furthermore, their models may provide inconsistent results if not every trial provides observations for the same time points. Lu et al. proposed a model for the cumulative incidence of a dichotomous measure over time but this model only allows events that can occur only once during the follow-up period [17]. In this paper, we present an alternative approach to network meta-analysis of repeated outcomes that will allow for an evaluation of the treatment effects over time that overcomes some of the limitations of the aforementioned models. The models we present assume a nonlinear development of treatment effects over time as described by fractional polynomials. Royston and Altman introduced fractional polynomials as an extension of polynomial models for determining the functional form of a continuous predictor [18]. The method is illustrated with an example based on RCTs evaluating interventions for osteoarthritis of the knee.

2. Network meta-analysis models for repeated measures using fractional polynomials If AB trials and AC trials are comparable in terms of effect modifiers (i.e. covariates that affect the treatment effect), then an indirect estimate for the treatment effect of C versus B (dBC ) can be obtained from the estimates of the effect of B versus A (dAB ) and the effect of C versus A (dAC ): dBC = dAC − dAB ; as such, transitivity holds [4]. In general, the treatment effect of intervention k relative to a control b can be expressed in terms of the overall reference treatment A ∶ dbk = dAk − dAb . When results are available at multiple time points for each treatment, a general random effects network meta-analysis model can be defined according to [9, 10]: { 𝜃jkt = 𝛿jbkt

𝜇jbt

if k = b, b ∈ {A,B,C}

𝜇jbt + 𝛿jbkt if k ≻ b ( ) ( ) ∼ normal dbkt , 𝜎t2 = normal dAkt − dAbt , 𝜎t2

Copyright © 2015 John Wiley & Sons, Ltd.

(1)

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

where 𝜃jkt is the linear predictor, a continuous measure such as ( log ) odds, change from baseline, or log hazard rate, for study j for treatment k at time point t. 𝜃jkt = g 𝛾jkt , where 𝛾jkt are the unknown parameters of the likelihood function for the data and g(.) is the link function, which maps the parameters of interest onto the plus/minus infinity range. 𝜇jbt are the study j specific outcomes at time t with comparator treatment b. 𝛿jbkt reflects the study and time-specific treatment effects of treatment k relative to comparator treatment b and is drawn from a normal distribution with the pooled estimates expressed in terms of the overall reference treatment A ∶ dbkt = dAkt − dAbt with dAAt = 0. Variance 𝜎t2 reflects the heterogeneity in treatment effects across studies for the different time points. If the heterogeneity is assumed to be constant over the different time points, the model as presented by Dakin et al. is obtained [16]. As an alternative to unconstrained baseline and treatment effects over time, we can make the assumption that outcomes can be described by a parametric function of time. A first-order fractional polynomial is obtained by describing the outcome of interest as a function of transformed time t in a linear model [18]: 𝜃t = 𝛽0 + 𝛽1 tp

(2)

The power p is chosen from the following set: −2. −1, −0.5, 0, 0.5, 1, 2, 3 with t0 = ln(t). Fractional polynomials are well suited for modeling nonlinear data and have been used in many applications, including survival analysis, meta-regression analysis, and meta-analysis of survival curves [19–26]. A second-order fractional polynomial is defined as: 𝜃t = 𝛽0 + 𝛽1 tp1 + 𝛽2 tp2 (3) If p1 = p2 = p, the model becomes a ’repeated powers’ model: 𝜃t = 𝛽0 + 𝛽1 tp + 𝛽2 tp ln t

(4)

Royston and Altman showed that by varying p1 and p2 and the parameters 𝛽0, 𝛽1 , and 𝛽2 , a wide range of curve shapes can be obtained [18, 19, 21, 23, 24]. For a network meta-analysis model with the outcome of treatment x as a parametric function of time, fx (t) where x = A, B, or C and t represents time, the consistency (assumption regarding the relative ) ( ) (t) − f (t) = f (t) − fA (t) − treatment effects underlying network meta-analysis translates into: f C B C ( ) fB (t) − fA (t) . Fractional polynomials can be incorporated in the network meta-analysis framework [25, 26]. The random effects model for a network meta-analysis of repeated outcomes over time based on a fractional polynomial of order M for k treatments labeled A, B, C, and so on can be described as:

𝜃jkt

⎧ M ⎪ 𝛽 + ∑ 𝛽 tpm with t0 = ln(t) if p1 ≠ .. ≠ pM 0jk mjk ⎪ m=1 =⎨ M ⎪ 𝛽 + 𝛽 tp1 + ∑ 𝛽 tp1 (ln(t))m−1 if M > 1, p = .. = p 0jk 1jk mjk 1 M ⎪ m=2 ⎩

⎧ ⎛ 𝜇0jb ⎞ ⎪⎜ ⋮ ⎟ if k = b, b ∈ {A,B,C} ⎟ ⎜ ⎛ 𝛽0jk ⎞ ⎪ ⎝ 𝜇Mjb ⎠ ⎜ ⋮ ⎟=⎪ ⎜𝛽 ⎟ ⎨ ⎝ Mjk ⎠ ⎪ ⎛ 𝜇0jb ⎞ ⎛ 𝛿0jbk ⎞ ⎪⎜ ⋮ ⎟+⎜ ⋮ ⎟ if k ≻ b ⎪⎜ ⎟ ⎜ ⎟ ⎩ ⎝ 𝜇Mjb ⎠ ⎝ 𝛿Mjbk ⎠ ⎛⎛ d0Ak ⎞ ⎛ d0Ab ⎞ ⎞ ⎛ 𝛿0jbk ⎞ ⎜ ⋮ ⎟ ∼ normal ⎜⎜ ⋮ ⎟ − ⎜ ⋮ ⎟ , Σ⎟ ⎜⎜ d ⎜𝛿 ⎟ ⎟ ⎜ ⎟ ⎟ ⎝⎝ MAk ⎠ ⎝ dMAb ⎠ ⎠ ⎝ Mjbk ⎠

(5)

2 · · · 𝜎0 𝜎M 𝜌0M ⎞ ⎛ 𝜎00 ⎟ ⋮ ⋱ ⋮ Σ=⎜ ⎜ ⎟ 2 𝜎 𝜎 𝜌 · · · 𝜎 ⎝ 0 M 0M ⎠ M

Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

where 𝜃jkt reflects the ‘underlying’ outcome in study j for treatment k at time point t and the link func⎛ 𝜇0jb ⎞ tion to transform this outcome to a normally distributed scale. The vectors ⎜ ⋮ ⎟ are trial-specific and ⎜𝜇 ⎟ ⎝ Mjb ⎠ ⎛ 𝛿0jbk ⎞ reflect the parameters 𝛽0 , 𝛽1 , … , 𝛽M of the comparator treatment, whereas the vectors ⎜ ⋮ ⎟ reflect ⎜𝛿 ⎟ ⎝ Mjbk ⎠ the study-specific difference in 𝛽0 , 𝛽1 , … , 𝛽M of the development of outcomes over time for treatment k relative to comparator treatment b and are drawn from a multivariate normal distribution with the pooled ⎛ d0AA ⎞ ⎛0⎞ estimates expressed in terms of the overall reference treatment A with ⎜ ⋮ ⎟ = ⎜ ⋮ ⎟. For example, ⎜d ⎟ ⎜0⎟ ⎝ MAA ⎠ ⎝ ⎠ ⎛ d0BC ⎞ ⎛ d0AC ⎞ ⎛ d0AB ⎞ ⎛ d0BD ⎞ ⎛ d0AD ⎞ ⎛ d0AB ⎞ ⎜ ⋮ ⎟ = ⎜ ⋮ ⎟ − ⎜ ⋮ ⎟, ⎜ ⋮ ⎟ = ⎜ ⋮ ⎟ − ⎜ ⋮ ⎟, and so on. Σ is the covariance matrix ⎜d ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ MBC ⎠ ⎝ dMAC ⎠ ⎝ dMAB ⎠ ⎝ dMBD ⎠ ⎝ dMAD ⎠ ⎝ dMAB ⎠ to reflect between-study heterogeneity regarding treatment effect parameters 𝛿mjbk that describe the development of treatment effects over time, which is assumed constant for all treatment comparisons. 𝜎m represent the variance for 𝛿mjbk (i.e. the difference in 𝛽m ) and 𝜌01 , 𝜌02 , … , 𝜌M−1,M is the correlation between these different treatment effect parameters of the polynomial. Of key interest from the analyses are the pooled estimates of dmAk and estimates for the heterogeneity. Please note that the relative treatment effects are changing over time once dm≥1 is different from 0. Under a fixed-effects model, the multivariate normal distribution with the pooled estimates will be ⎛ 𝛿0jbk ⎞ ⎛ d0Ak ⎞ ⎛ d0Ab ⎞ replaced with ⎜ ⋮ ⎟ = ⎜ ⋮ ⎟ − ⎜ ⋮ ⎟ and as a result, the between-study covariance matrix does ⎜𝛿 ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ Mjbk ⎠ ⎝ dMAk ⎠ ⎝ dMAb ⎠ not need to be estimated. When heterogeneity is assumed only for d0Ak and the other effect parame⎛⎛ d0Ak ⎞ ⎛ d0Ab ⎞ ⎞ ⎛ 𝛿0jbk ⎞ ters d1Ak , … , dMAk are fixed, then ⎜ ⋮ ⎟ ∼ normal ⎜⎜ ⋮ ⎟ − ⎜ ⋮ ⎟ , Σ⎟ is replaced with 𝛿0jbk ∼ ⎜⎜ d ⎜𝛿 ⎟ ⎟ ⎜ ⎟ ⎟ ⎝⎝ MAk ⎠ ⎝ dMAb ⎠ ⎠ ⎝ Mjbk ⎠ ⎛ 𝛿1jbk ⎞ ⎛ d1Ak ⎞ ⎛ d1Ab ⎞ ( ) normal d0Ak − d0Ab , 𝜎 2 and ⎜ ⋮ ⎟ = ⎜ ⋮ ⎟ − ⎜ ⋮ ⎟. ⎜𝛿 ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ Mjbk ⎠ ⎝ dMAk ⎠ ⎝ dMAb ⎠ A random effects model with only a heterogeneity parameter for d0Ak implies that the between-study variance of effect estimates remains constant over time. Random effects models with (additional) heterogeneity parameters for d1Ak , … , dMAk have the flexibility to capture between-study variance regarding changes in the relative treatment effects over time. The random effects fractional polynomial model presented earlier does not account for correlation between trial-specific 𝛿s in multiple-arm trials (>2 treatments). Bayesian random effects fractional polynomials models with only a heterogeneity parameter for d0Ak can be easily extended to fit trials with three or more treatment arms by decomposition of a multivariate normal distribution ⎛ 𝛿0jbk1 ⎞ as a series of conditional univariate distributions, as shown by Cooper et al. [27]. If ⎜ ⋮ ⎟ ∼ ⎜𝛿 ⎟ ⎝ 0jbkp ⎠ 2 𝜎 2∕2 · · · 𝜎 2∕2 𝜎 ⎛⎛ d ⎞⎞ ⎞ ⎛ ⎜⎜ 0bk1 ⎟ ⎜ 𝜎 2∕2 𝜎 2 · · · 𝜎 2∕2 ⎟⎟ normal ⎜ ⋮ , ⎜ , then the conditional univariate distributions for arm i given the ⋮ ⋮ ⋱ ⋮ ⎟⎟ ⎜⎜⎝ d0bkp ⎟⎠ ⎜ 𝜎 2 𝜎 2 ⎟⎟ 2 ⎝ ⎝ ∕2 ∕2 · · · 𝜎 ⎠⎠ previous 1,….(i-1) arms are ) ( |⎛ 𝛿 i−1 ( ) (i + 1) | 0jbk1 ⎞ ∑ 1 |⎜ 2 ⋮ ⎟ ∼ normal d0bki + 𝛿0jbki | 𝜎 − d0bkj , 𝛿 |⎜ i j=1 0jbkj 2i |⎝ 𝛿0jbki−1 ⎟⎠ | Different values for the powers pm of the fractional polynomials correspond to different models. The best fitting model can be selected based on clinical expertise regarding the expected development Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

of treatment effects over time and goodness-of-fit comparisons. The goodness-of-fit can be computed as the difference between the deviance for the fitted model and the deviance for the saturated model (which fits the data perfectly). Within a frequentist framework, the Akaike information criterion can be used for model selection [28]. In a Bayesian framework, the Bayesian information criterion or deviance information criterion (DIC) can be used [29, 30]. Model 5 is a general representation of the fractional polynomial network meta-analysis model for repeated measures. Depending on the evidence base, simplifications may be required to ensure sufficient data points are available to obtain stable model parameter estimates. Typically, a second-order fractional polynomial is sufficient to capture changes in outcomes over time [22]. Further simplifications can be made when the outcome of interest reflects change from baseline by removing the intercept and associated treatment-effect parameters (d0Ak )

3. Illustrative example An example of the repeated measures models is presented for a network meta-analysis of treatment for osteoarthritis (OA) of the knee with different hyaluronan (HA)-based viscosupplements. OA is defined as a painful degenerative process that is not primarily inflammatory, involving progressive deterioration of articular cartilage as well as reactive alterations of subchondral bone and the surrounding joint structure, which includes the synovial lining, the ligaments, the tendons, the capsule, and the periarticular bone [31]. Symptomatic knee OA results in joint pain, dysfunction, and stiffness related to inactivity [32,33]. Viscosupplementation with HA-based products is a therapy for the treatment of OA, based on the importance of HA in synovial fluids. Injection of HA-based viscosupplements in the joint helps restore visco-elasticity of the HA, thereby decreasing pain, among other symptoms [34–38]. Sixteen RCTs investigating analgesia with HA products were identified with a systematic literature search [39–54]. The outcome measure of interest was change from baseline (CFB) in pain at movement (0 − 100 mm visual analogue scale). The different RCTs evaluated the following interventions: three, four, or five injections with HA with a molecular weight (MV) of 0.5 − 0.73 million Da (Hyalgan) (3 Hy−0.5−0.73; 4 Hy−0.5−0.73; 5 Hy−0.5−0.73); three injections with HA MW of 0.62−11.7 million Da (Supartz) (3 Hy−0.62 − 11.7); three injections with HA MW of 2.4 − 3.6 million Da (Euflexxa) (3 Hy−2.4 − 3.6); and three injections of Hylan GF-20 MW 6 million Da (Synvisc) (3 HyGF20). The identified trials are presented in the evidence network in Figure 1. The CFB in pain, along with standard errors for each intervention over time as obtained from the RCTs, is presented in Table I. The CFB in pain at each time point are correlated over time. Unfortunately, this within-trial correlation was not reported for the included studies. As such, we performed sensitivity analyses assuming different values for the correlation: 0, 0.5, and 0.9. In addition, an analysis was performed assuming a uniform distribution ranging from 0 to 0.95 for the correlation. In order to properly incorporate the within-trial correlation by treatment group, a multivariate normal likelihood is required. Because the different studies have different numbers of follow-up measures, we decided to acknowledge within-trial correlation with a more practical approach by using a univariate normal likelihood distribution

Figure 1. Network of randomized controlled trials.

Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

Week 4

Week 5

Copyright © 2015 John Wiley & Sons, Ltd.

3 HY−2.4 − 3.6 3 HYGF20

3 HY−0.5 − 0.73 5 HY−0.5 − 0.73

Stitik et al. (2007)

3 HYGF20 PLACEBO

Wobig et al. (1998)

Kirchner et al. (2006)

3 HYGF20 PLACEBO

Scale et al. (1994)

3 HY−0.62 − 11.7 3 HYGF20 PLACEBO

Karlsson et al. (2002)

−38.3 −21.3

Huskisson & Donnelly (1999) 5 HY−0.5 − 0.73 PLACEBO

3 HY−0.5 − 0.73 PLACEBO

−15.6 −8.7 −14.2 −18.0

−0.9 −5.6

−12.0 −5.0

−9.5 −2.5

−5.0 −7.0 −7.0

5.2 7.4

−5.7 −15.2

−26.0 −16.0

−32.0 −21.0

5.8 4.2 4.6 4.2

−12.0 −16.0 −11.0

2.0 2.0 3.2

5.9 8.0

4.6 4.2

5.8 4.2

2.6 2.5 3.6

−39.0 −22.0

−44.0 −20.0

−20.0 −18.0 −21.0

−31.0 −25.0

4.6 4.2

5.0 4.2

2.8 2.9 4.0

4.0 3.0

2.5 2.8

5 HY−0.5 − 0.73 5 HY−0.5 − 0.73 PLACEBO PLACEBO

3.0 2.0

2.4 3.7

−24.3 −16.2

Henderson et al. (1994)

−20.0 −15.0

−18.3 −13.0

−12.5 −11.0

−33.0 −26.0

1.9 1.9

5.0 3.0

−13.2 −20.5

−33.3 −20.7

−21.5 −18.1 −8.0

−29.0 −30.0

4 HY−0.5 − 0.73 PLACEBO

2.0 3.9

4.2 4.1 3.3

2.9 3.1

Dougados et al. (1993)

−6.3 −9.0

−17.3 −14.6 −8.4

−27.0 −27.0

3 HYGF20 PLACEBO

2.9 3.1

Cubukcu et al. (2005)

−23.0 −21.0

5 HY−0.5 − 0.73 PLACEBO

4.0 4.3 3.3

3.0 3.1

Corrado et al. (1995)

−8.3 −6.0 −7.2

−19.0 −20.0

5 HY−0.5 − 0.73 3 HY−0.5 − 0.73 PLACEBO

3.1 3.1

Carrabba et al. (1995)

Jubb et al. (2003)

Week 7

Week 8

Week 9

Week 11

Week 12

5.9 8.0

4.1 5.3

2.6 1.9 2.2 1.7

7.5 7.2

4.2 4.6 3.6

2.9 3.0

−35.5 −25.8

3.6 2.9

−46.0 −15.0

−52.0 −19.0

−33.5 −19.8

−31.0 −13.2

−39.0 −19.1

−21.2 −16.3 −4.6

4.6 4.2

5.8 4.2

4.5 5.3

2.4 3.1

6.5 6.7

3.9 4.9 3.8

−30.0 −32.0

6.1 6.1

−10.0 −5.0

2.1 1.9

−23.2 −24.3

−31.2 −28.7

−46.0 −12.0

−53.0 −17.0

−22.0 −22.0 −19.0

−31.0 −31.0

5.9 8.0

2.0 2.0

4.6 4.2

3.6 3.0

3.2 3.5 4.6

3.0 3.0

CFB SEM CFB SEM CFB SEM CFB SEM CFB SEM CFB SEM CFB SEM CFB SEM CFB SEM CFB SEM

Week 3

5 HY−0.5 − 0.73 −15.0 PLACEBO −15.0

Intervention

Week 2

Altman et al. (1998)

Study

Week 1

Table I. Individual study results; change from baseline regarding pain on movement.

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Statist. Med. 2015

Copyright © 2015 John Wiley & Sons, Ltd.

3 HYGF20 PLACEBO

5 HY−0.5 − 0.73 PLACEBO

5 HY−0.5 − 0.73 3 HY−0.5 − 0.73 PLACEBO

5 HY−0.5 − 0.73 PLACEBO

3 HYGF20 PLACEBO

4 HY−0.5 − 0.73 PLACEBO

5 HY−0.5 − 0.73 5 HY-0.5-0.73 PLACEBO PLACEBO

5 HY−0.5 − 0.73 PLACEBO

3 HY−0.5 − 0.73 PLACEBO

Dickson et al. (2001)

Altman et al. (1998)

Carrabba et al. (1995)

Corrado et al. (1995)

Cubukcu et al. (2005)

Dougados et al. (1993)

Henderson et al. (1994)

Huskisson & Donnelly (1999)

Jubb et al. (2003)

3 HY−0.62 − 11.7 3 HYGF20 PLACEBO

3 HY−2.4 − 3.6 PLACEBO

Altman et al. (2009)

Karlsson et al. (2002)

3 HYGF20 3 HY−0.62 − 11.7

Intervention

Wobig et al. (1999)

Study

Table I. Continued.

−32.8 −13.6

−33.0 −33.0

−13.5 −10.0

−8.0 −10.0

CFB

4.9 5.5

2.9 3.0

1.0 1.0

3.6 4.5

SEM

Week 16

−9.0 −10.0

−20.0 −17.5

−18.0 −20.0

CFB

2.2 1.9

1.0 1.0

3.6 4.5

SEM

Week 18

−21.0 −27.0 −19.0

−24.0 −22.5

−29.0 −30.0

CFB

3.2 3.5 4.2

1.0 1.0

4.2 4.5

SEM

Week 20

−37.0 −34.0

CFB

2.8 2.9

SEM

Week 21

−16.0 −20.0 −21.0

−26.4 −8.2

−36.0 −31.0

CFB

3.8 3.7 4.5

4.7 5.3

2.8 3.1

SEM

Week 26

−25.0 −22.0

CFB

1.0 1.0

SEM

Week 36

−38.9 −32.7

−37.0 −30.0

CFB

4.2 3.9

4.2 4.5

SEM

Week 52

−39.0 −26.0

−25.5 −22.0

−37.0 −20.0

4.0 4.0

1.0 1.0

4.2 4.5

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Statist. Med. 2015

Copyright © 2015 John Wiley & Sons, Ltd.

3 HYGF20 PLACEBO

3 HYGF20 PLACEBO

3 HY−2.4 − 3.6 3 HYGF20

3 HY−0.5 − 0.73 5 HY−0.5 − 0.73

3 HYGF20 3 HY−0.62 − 11.7

3 HY−2.4 − 3.6 PLACEBO

3 HYGF20 PLACEBO

Scale et al. (1994)

Wobig et al. (1998)

Kirchner et al. (2006)

Stitik et al. (2007)

Wobig et al. (1999)

Altman et al. (2009)

Dickson et al. (2001)

CFB

SEM

Week 16

*CFB = change from baseline; *SEM = standard error of the mean

Intervention

Study

Table I. Continued.

−27.5 −21.0

CFB

1.0 1.0

SEM

Week 18 CFB

SEM

Week 20 CFB

SEM

Week 21

−25.0 −17.5

−21.8 −23.3

CFB

1.0 1.0

5.9 8.0

SEM

Week 26

−10.4 −26.3

CFB

5.9 8.0

SEM

Week 36

−16.1 −35.4

CFB

5.9 8.0

SEM

Week 52

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

for the CFB at every time point with adjustment of the variance: 2 cfbjkt ∼ normal(𝜃jkt , 𝜎jkt )

(6)

2 where cfbjkt is the observed CFB in pain at time point t with treatment k in study j. 𝜎jkt is the correspond2 ing variance as a measure of uncertainty. Estimates for 𝜎jkt were obtained from the reported standard errors as presented in Table I and adjusted by dividing these variance estimates by 1 − 𝜌2 , where 𝜌 is the assumed correlation coefficient between subsequent time points. With this approach, the correlation between subsequent time points is assumed constant over time. 𝜃jkt is the true underlying CFB in pain at time point t with treatment k in study j. Given the general pattern in outcomes observed in the individual studies, it was decided that a second-order fractional polynomial network meta-analysis model was sufficiently flexible to capture the development of treatment effects over time. The following model was used for the evidence synthesis: { p1 ≠ p2 𝛽1jk tp1 + 𝛽2jk tp2 with t0 = ln(t) 𝜃jkt = 𝛽1jk tp1 + 𝛽2jk tp1 ln(t) p1 = p2

(

𝛽1jk 𝛽2jk

)

𝛿1jbk (

d1AA d2AA

) ( ⎧ 𝜇1jb if k = b, b ∈ {A, B, C, D, E, F} ⎪ 𝜇 ) = ⎨ ( 2jb ) ( 𝜇1jb 𝛿1jbk ⎪ + if k ≻ b ⎩ 𝜇2jb d2Ak − d2Ab ( ) ( ) ∼ normal d1bk , 𝜎12 = normal d1Ak − d1Ab , 𝜎12

(7)

) =0

where 𝜃jkt , the true underlying CFB in pain at time point t with treatment k in study j , is now described as a function of time t with p = {−2, −1, −0.5, 0, 0.5, 1, 2, 3} and t0 = ln(t) with treatment and studyspecific shape parameters 𝛽1jk and 𝛽2jk . Because the endpoint of interest is CFB, 𝛽ojk from equation 5 equals 0 and therefore drops out of the model. In the random effects models used for the network metaanalysis in this example, there is one between-study heterogeneity parameter, 𝜎12 , related to the relative treatment effects for 𝛽1jk . As an alternative, the between-study heterogeneity was assumed to concern only treatment effects in terms of 𝛽2jk in a separate analysis. Fixed effects models were also evaluated. Here, the prior distributions used for the meta-analysis are presented. ( ) (( ) ) ) ( 4 𝜇1jb 0 10 0 ∼ normal , T𝜇 T𝜇 = 0 𝜇2jb 0 104 ( ) (( ) ) ) ( 4 d1Ak 0 10 0 ∼ normal , Td Td = 0 d2Ak 0 104 𝜎1 ∼ uniform(0, 10) In the alternative analysis with study heterogeneity assumed to concern only treatment effects in terms of 𝛽2jk , the prior 𝜎1 ∼ uniform(0, 10) was replaced with 𝜎2 ∼ uniform(0, 10). With the fixed effects model, it is not necessary to define prior distributions for 𝜎1 or 𝜎2 . The parameters of the different models were estimated using a Markov Chain Monte Carlo method, as implemented in the OpenBUGS software package [55] (See Appendix for the code). The first 20,000 iterations from the sampler were discarded as ‘burn-in’, and the inferences were based on an additional 80,000 iterations using two chains. Convergence of the chains was confirmed by the Gelman– Rubin statistic. The DIC was used to compare the goodness-of-fit of the different models [29, 30]. DIC provides a ̂ [24]. measure of model fit that penalizes model complexity according to DIC = D + pD, pD = D − D ̂ is the deviance D is the posterior mean deviance [30], pD is the ‘effective number of parameters’, and D evaluated at the posterior mean of the model parameters. In general, a more complex model will result in a better fit to the data, demonstrating a smaller residual deviance. The model with the lowest DIC is Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Table II. Goodness-of-fit estimates for fixed effects and random effects fractional polynomial network metaanalysis models with different powers p1 and p2 .

Fixed effects model p1

−2 −2 −2 −2 −2 −2 −2 −2 −1 −1 −1 −1 −1 −1 −1 −0.5 −0.5 −0.5 −0.5 −0.5 −0.5 0 0 0 0 0 0.5 0.5 0.5 0.5 1 1 1 2 2 3

Random effects model with heterogeneity for treatment effect for Beta1

Random effects model with heterogeneity for treatment effect for Beta2

p2

Dbar

pD

DIC

Dbar

pD

DIC

Dbar

pD

DIC

−2 −1 −0.5 0 0.5 1 2 3 −1 -0.5 0 0.5 1 2 3 −0.5 0 0.5 1 2 3 0 0.5 1 2 3 0.5 1 2 3 1 2 3 2 3 3

7193 4735 2477 1855 1984 3730 5916 6978 2828 1509 1459 1708 2842 4465 5343 1039 1222 1463 2057 2958 3479 1374 2128 1831 1798 1931 1037 1017 1170 1321 1278 2072 2566 3692 4467 5327

36.2 38.7 39.9 40.5 40.7 40.7 40.8 40.8 40.7 40.4 41.6 41.7 41.7 41.9 41.8 41.4 41.9 41.8 42.0 41.9 42.0 41.9 41.0 42.0 41.9 42.0 42.0 41.9 41.9 42.0 42.0 42.1 41.9 42.0 42.1 42.1

7229 4773 2517 1895 2024 3770 5956 7019 2868 1549 1501 1750 2884 4506 5385 1081 1263 1504 2099 3000 3521 1416 2169 1873 1840 1973 1079 1059 1211 1363 1320 2114 2608 3734 4510 5369

7188 4724 2460 1829 1958 3712 5909 6975 2790 1463 1407 1657 2799 4433 5314 944 1123 1367 1969 2882 3406 1229 1965 1687 1653 1783 886 874 1031 1181 1145 1935 2422 3567 4337 5210

38.2 40.6 41.7 42.8 42.7 43.2 43.4 43.5 45.9 45.9 47.4 47.4 47.4 47.7 47.7 50.1 50.5 50.5 50.6 50.8 50.6 53.5 52.9 53.4 53.6 53.6 53.7 53.9 53.6 53.6 53.6 53.5 53.3 52.3 53.3 52.0

7226 4765 2502 1871 2001 3755 5952 7019 2836 1509 1454 1704 2846 4480 5362 995 1173 1417 2019 2933 3457 1282 2018 1741 1707 1836 940 928 1085 1234 1198 1988 2476 3619 4391 5262

7185 4703 2396 1692 1819 3557 5744 6816 2736 1420 1298 1543 2672 4297 5186 893.9 1061 1300 1892 2801 3334 1236 1966 1693 1677 1824 891.3 877.9 1051 1218 1153 1953 2461 3577 4359 5220

37.0 44.0 48.4 52.5 52.5 52.5 52.4 52.1 46.4 48.9 53.6 53.4 53.6 53.5 53.2 51.8 53.5 53.7 53.7 53.5 53.3 53.5 52.7 53.4 53.3 52.9 53.5 53.5 53.5 53.0 53.2 53.5 52.8 53.1 53.2 51.2

7222 4747 2444 1744 1872 3609 5796 6868 2783 1469 1351 1597 2726 4351 5239 946 1115 1354 1946 2855 3388 1289 2019 1747 1730 1877 945 932 1105 1271 1206 2007 2514 3630 4412 5272

Dbar, posterior mean deviance; pD, the effective number of parameters as a measure of model complexity; DIC, Deviance Information Criteria

the model providing the ‘best’ fit to the data adjusted for the number of parameters. In Table II, the fit of the different models to the data (16 RCTs, 36 intervention arms, and 152 data points) are summarized with model fit statistics. The random effects models were associated with a smaller residual deviance than the fixed effects models. Taking into account the increased model complexity of the random effects approach, the DIC with the random effects models was also lower. Hence, we prefer to use the random effects models over the fixed effects model. Of all the random effects models, those with p1 = 0.5 and p2 = 1 were associated with the smallest DIC. The model with a heterogeneity parameter for d1Ak was slightly more favorable ( (∼3)points) than the model with a heterogeneity parameter for d2Ak . d1bk corresponding to the random effects model of choice (p1 = 0.5, p2 = 1) are The estimates for d2bk presented in Tables III (d1bk ) and IV (d2bk ) below the diagonal. These relative treatment effects correspond to interventions in each row, relative to the intervention in the column. Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Table III. Model parameter estimates for d1Ak as obtained with random effects fractional polynomial network meta-analysis model with powers p1 = 0.5 and p2 = 1 (white cells) and corresponding independent means model for direct effects (grey cells). Estimates are expressed as median (and 95%credible interval) of posterior distributions.

Estimates presented in white cells obtained from network meta-analysis model are relative treatment effects of interventions in rows relative to comparator in columns. Estimates presented in grey cells obtained from independent means model are relative treatment effects of interventions in columns relative to comparator in rows.

Table IV. Model parameter estimates for d2Ak as obtained with random effects fractional polynomial network meta-analysis model with powers p1 = 0.5 and p2 = 1 (white cells) and corresponding independent means model for direct effects (grey cells). Estimates are expressed as median (and 95%credible interval) of posterior distributions.

Estimates presented in white cells obtained from network meta-analysis model are relative treatment effects of interventions in rows relative to comparator in columns. Estimates presented in grey cells obtained from independent means model are relative treatment effects of interventions in columns relative to comparator in rows.

With a network meta-analysis, the assumption is made that relative treatment effects between interventions based on direct comparisons are consistent with estimates of relative treatment effects based on indirect comparisons for the same contrasts. As mentioned before, inconsistency occurs if there are systematic differences in relative treatment effect modifiers between different types of direct comparisons. In this example, there is both direct and indirect evidence for some treatment comparisons (Figure 1), which allows for a comparison of the treatment effect estimates based on direct evidence only with those obtained from the network meta-analysis model. Here, an independent means model according to Dias et al. is presented that synthesizes only direct evidence without making the assumption that direct and indirect evidence are consistent [56]. Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

{ 𝜃jkt =

𝛽1jk tp1 + 𝛽2jk tp2 𝛽1jk

tp1

+ 𝛽2jk

tp1

p1 ≠ p2 ln(t)

p1 = p2

with t0 = ln(t)

⎧(𝜇 ) if k = b, b ∈ {A, B, C, D, E, F} ( ) ⎪ 𝜇1jb ⎪ 2jb 𝛽1jk = ⎨( ) ( ) 𝛽2jk ⎪ 𝜇1jb + 𝛿1jbk if k ≻ b ⎪ 𝜇2jb d2bk ⎩ ( ) 𝛿1jbk ∼ normal d1bk , 𝜎12

(8)

In this model, each of the contrasts for which direct evidence is available represents a separate, unrelated, basic parameter to be estimated. The network meta-analysis (equation 7) estimates six basic parameters for relative treatment effects corresponding to 𝛽1 (d1AB , d1AC , d1AD , d1AE , d1AF , and d1AG ) and 𝛽2 (d2AB , d2AC , d2AD , d2AE , d2AF , and d2AG ), whereas the independent-means model estimates nine unrelated relative treatment effect parameters 𝛽1 (d1AB , d1AC , d1AD , d1AE , d1AF , d1AG , d1BD , d1EG , and d1FG ) and 𝛽2 (d2AB , d2AC , d2AD , d2AE , d2AF , d2AG , d2BD , d2EG , and d2FG ). The prior distributions for the independent means model are presented here: ) (( ) ) ) ( ( 4 𝜇1jb 0 10 0 ∼ normal , T𝜇 T𝜇 = 0 𝜇2jb 0 104 ( ) (( ) ) d1bk 0 ∼ normal , Td 0 d2bk ) ( 4 10 0 𝜎1 ∼ uniform(0, 10) Td = 0 104 ) ( d1bk are presented above the diagonal in Tables III and IV. These indeThe parameter estimates for d2bk pendent relative treatment effects correspond to interventions in each column, relative to the comparator in the row. The estimates obtained with the network meta-analysis model (estimates below the diagonal) can be compared with the corresponding direct estimates above the diagonal. The comparison of 3 HyGF20 versus 3 Hy-2.4-3.6 shows an obvious difference between the estimated direct evidence from the independent means model in comparison with the estimated mixed effect from the network metaanalysis model. However, the only study that provides direct evidence for this treatment comparison is by Kirchner et al., which includes data at 12 weeks follow-up only [50]. As such, there is no sufficient data to estimate a curve of treatment effects over time based on direct evidence; we have to be very careful when interpreting the observed difference between the results for this treatment contrast obtained with

Figure 2. Treatment effects of the different interventions relative to placebo as obtained with the random effects fractional polynomial network meta-analysis model with powers p1 = 0.5 and p2 = 1, assuming no within-trial correlation of outcomes over time. The thick lines reflect the point estimates and the thin lines in the same color reflect the 95% credible intervals. Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Figure 3. Rankogram reflecting the probability that each intervention ranks first, second, third, and so on out of the seven interventions compared regarding the area under the pooled change from baseline curves up to 26 weeks of follow-up as obtained with the network meta-analysis.

the two models. Differences between estimates from the independent means model and network metaanalysis model were also observed for the comparisons between 3HyGF20 and 3 Hy-0.62-11.7, as well as between 5 Hy-0.5-0.73 and 3 Hy-0.5-0.73. However, given the large degree of uncertainty in the direct estimates for d1bk and d2bk with these comparisons, it is difficult to judge whether the observed differences in the results obtained with the independent means and network meta-analysis models are due to chance or the result of inconsistency between direct and indirect evidence combined in the mixed estimates of the network meta-analysis model. Given the treatment effect parameters obtained with the network meta-analysis model as presented in Tables III and IV, the corresponding treatment effects over time relative to placebo are presented in Figure 2. Based on the area under each curve over the 26-week period, the probability that each treatment was best, ranked first, second, third, fourth, fifth, sixth, and last was calculated and presented in Figure 3. The estimates for the basic parameters reflecting the treatment effects of each intervention relative to placebo (d1Ak and d2Ak ) for the different sensitivity analyses performed are presented in Table V. As expected, the stronger the assumed correlation, the greater the uncertainty in the parameter estimates. In addition, the shape of the curves changed quite a bit when a correlation of 0.9 was assumed instead of 0. In Figure 4, the corresponding treatment effects over time relative to placebo are presented. The sensitivity analysis assuming a uniform distribution for the correlation coefficient resulted in similar estimates as the analysis with a correlation of 0.5, but with slightly wider credible intervals reflecting the uncertainty in the assumed correlation.

4. Discussion In this paper, a method for network meta-analysis of repeated measures using a multi-dimensional treatment effect is presented. With fractional polynomials, the development of outcomes over time for the interventions compared in a trial is modeled, and the differences in the parameters of these fractional polynomials within a trial are considered the multidimensional treatment effect, which are synthesized (and indirectly compared) across studies. In contrast to previous methods for repeated measures, these evidence synthesis models include the evaluation of nonlinear trends in treatment effects and simultaneously analyze studies with different follow-up periods and outcomes measured at different time points. With multiple effect sizes from the same study correlated, the assumption of independence is violated, and results of a meta-analysis of repeated measures ignoring correlation may be distorted. When IPD is available for all trials, the correlation between repeated measures can be taken into account. Jones et al. have shown that ignoring correlation of outcomes over time in an AD meta-analysis has a different impact on the results, depending on whether time is considered a factor in the meta-analysis model or as a continuous variable [12]. With the first approach, ignoring the within-trial correlation results in overestimation of the uncertainty of pooled treatment effects at the different time points, whereas with the second approach, the uncertainty is underestimated. Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

Copyright © 2015 John Wiley & Sons, Ltd.

−2.14 −4.58 −1.68 −2.79 −3.59 −4.51 0.34 0.51 −0.07 0.08 −0.03 −0.28

3 Hy−0.5 − 0.73 (k =B) 4 Hy−0.5 − 0.73 (k =C) 5 Hy−0.5 − 0.73 (k =D) 3 Hy−0.62 − 11.7 (k =E) 3 Hy−2.4 − 3.6 (k =F) 3 HyGF20 (k =G)

3 Hy−0.5 − 0.73 (k =B) 4 Hy−0.5 − 0.73 (k =C) 5 Hy−0.5 − 0.73 (k =D) 3 Hy−0.62 − 11.7 (k =E) 3 Hy−2.4 − 3.6 (k =F) 3 HyGF20 (k =G)

(−0.54, 1.2) (−0.03, 1.05) (−0.70, 0.56) (−0.94, 1.10) (−0.33, 0.26) (−1.20, 0.65)

(−7.69, 3.43) (−12.97, 3.85) (−5.52, 2.02) (−9.41, 3.99) (−9.48, 2.25) (−8.93, 0.02)

𝝆=0

0.33 0.51 −0.10 0.10 −0.02 −0.25

−2.11 −4.53 −1.53 −2.84 −3.58 −4.59

d1Ak

(−0.65, 1.32) (−0.12, 1.13) (−0.83, 0.65) (−1.08, 1.30) (−0.38, 0.31) (−1.37, 0.82)

0.36 0.51 −0.17 0.51 −0.02 0.34

(−7.90, 3.64) −2.04 (−12.98, 3.91) −4.54 (−5.60, 2.35) −0.94 (−9.86, 4.21) −3.99 (−9.47, 2.22) −2.85 (−9.33, 0.33) −5.99 d2Ak

𝝆 = 0.5

(−1.44, 2.14) (−0.76, 1.77) (−1.50, 1.19) (−1.83, 2.71) (−0.71, 0.64) (−1.71, 2.38)

(−9.81, 5.82) (−14.28, 5.30) (−6.49, 4.40) (−13.43, 5.91) (−8.69, 2.38) (−13.08, 1.29)

𝝆 = 0.9

0.36 0.51 −0.08 0.13 −0.03 −0.21

−2.17 −4.55 −1.53 −2.96 −3.51 −4.69

(−0.76, 1.50) (−0.20, 1.22) (−0.90, 0.75) (−1.17, 1.46) (−0.42, 0.37) (−1.39, 1.00)

(−8.25, 3.89) (−13.17, 4.02) (−5.75, 2.55) (−10.20, 4.47) (−9.38, 2.27) (−9.71, 0.46)

𝝆 ∼ 𝐮𝐧𝐢𝐟 𝐨𝐫𝐦(0, 0.95)

Table V. Sensitivity analyses regarding the impact of within-trial correlation (𝜌) on model parameter estimates for d1Ak and d2Ak as obtained with random effects fractional polynomial network meta-analysis model with powers p1 = 0.5 and p2 = 1. Estimates are expressed as median (and 95%credible interval) of posterior distributions.

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

Figure 4. Treatment effects of the different interventions relative to placebo as obtained with the random effects fractional polynomial network meta-analysis model with powers p1 = 0.5 and p2 = 1 assuming a within-trial correlation of 0.9 for outcomes at subsequent time points. The thick lines reflect the point estimates and the thin lines in the same color reflect the 95% credible intervals.

Because publications typically do not report the correlation or covariance matrices of outcomes over time, meta-analysis based on AD requires assumptions about the correlations [12–14]. Given the continuous nature of the outcome of interest in our example, a multivariate normal likelihood would be appropriate to capture the assumed correlation over time. However, we opted for a pragmatic approximation by using univariate normal likelihood distributions, which effectively ignored covariance but inflated the variance at each time point to match the variation, as would have been observed in a multivariate normal likelihood distribution. This approach implicitly assumes that the correlation between subsequent time points remains constant over time. Nonetheless, the use of multivariate likelihoods is a more elegant approach, which would allow for the evaluation of different covariance structures in relation to time. Perhaps the best general approach to estimate within-trial correlation of repeated measures over time is the use of patient level data. However, patient level data are often not available for all relevant studies in the network. If for select studies within a network, patient level data are available, or to report withinstudy correlations, a multivariate normal distribution could be used for those studies, and a univariate normal distribution could be used for the remaining studies, where the parameters related to the treatment differences over time and the between-study heterogeneity would be shared. An easier approach is to obtain the within-study correlation from the patient level data set or external evidence and impute this in the study level model for the trials with missing information [57]. Similar network meta-analysis models have been presented for the synthesis of published time-toevent or survival curves [25, 26]. However, analyses of these outcomes do not require assumptions about autocorrelations if the meta-analysis is based on the underlying hazard function over time. The hazards at a certain time point are conditional upon not having experienced the event up to that time point and therefore are uncorrelated. As an alternative to fractional polynomials, spline functions can also be used to describe outcomes over time. Although a wide range of curve shapes can be obtained with (second order) fractional polynomials, splines can be more flexible, depending upon the number of nodes [18–24]. However, for the consistency assumption to hold in the network meta-analysis framework, the nodes of (cubic) splines need to be set at the same locations across all arms of all trials, which will restrict the possible shapes, flexibility, and possibly the model fit. Based on our experience, network meta-analysis of repeated measures with fractional polynomials provides a more practical compromise between flexibility and efficiency regarding model estimation than network meta-analysis with cubic splines. With the models introduced in this paper, we need to ensure there are a sufficient number of time points available within studies to ensure stable parameter estimates. To minimize the risk of treatment effect patterns over time that are affected by bias similar to ecological inference fallacy, one needs for each intervention in the network at least one study-arm with more data points than parameters describing the development of the outcome over time. However, additional research is needed to better understand the interaction between available time points, the available direct and indirect treatment comparisons at the different time points, model complexity, and the associated risk of spurious or confounded findings. A simulation study can be an excellent approach to better understand the performance of these models, Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE

given different assumptions regarding the data structure. One approach to select relevant parametric models, whether fractional polynomials or splines, is to start with a nonparametric saturated model, as done by Dakin et al., and to compare the more structured models with the saturated model [16]. However, it is important to consider the clinical rationale and to define potentially relevant expected patterns of treatment effects over time prior to the analyses in order to guide model selection. This is especially important when substantial between-study heterogeneity is present and many different competing models may be equally valid from a statistical perspective. Apart from RCTs, in many disease areas, the evidence base also includes noncontrolled trials or single group observational evidence. These studies can be very useful to incorporate in the evidence synthesis of repeated measures with fractional polynomials, not to estimate relative treatment effects between interventions but to improve estimation of the (nonlinear) trend of outcomes over time. This is especially the case if the uncontrolled studies have a relatively large sample size and a longer follow-up than the RCTs. A limitation of the fractional polynomial models is the interpretation of the treatment effect parameters. We provided these estimates in Tables III and IV to evaluate consistency between the network meta-analysis and independent means models, but for interpretation of results, it is more useful and recommended to plot treatment effects over time (Figure 2) or present these at multiple time points in a table. The Bayesian approach to the analysis allows ranking of treatments based on the location, spread, and overlap of the posterior distributions of the treatment effect parameters. Given that multiple parameters describe the development of treatment effects over time for each intervention, ranking treatments according to each of these parameters is not informative. In the example, we opted to use the area under the treatment effect curves to obtain one measure used for the ranking of treatments. The area under the curve reflects the average of the treatment effects over time. Alternatively, one can also rank the treatments at specific time points, but this may complicate decision-making because multiple ‘rankings’ will be obtained. In conclusion, network meta-analysis with models that represent the treatment effects in terms of several parameters using fractional polynomials can be considered a useful addition to models for network meta-analysis of repeated measures previously proposed. When the outcomes of interest across different trials are not reported at the same time points, the proposed models provide clear advantage over existing models.

References 1. Inthout J, Ioannidis JP, Borm GF. Obtaining evidence by a single well-powered trial or several modestly powered trials. Statistical Methods in Medical Research 2012. DOI: 10.1177/0962280212461098. 2. Higgins JPT, Whitehead A. Borrowing strength from external trials in a meta-analysis. Statistics in Medicine 1996; 15:2733–2749. 3. Ades AE. A chain of evidence with mixed comparisons: models for multi-parameter synthesis and consistency of evidence. Statistics in Medicine 2003; 22(19):2995–3016. 4. Lu G, Ades AE. Combination of direct and indirect evidence in mixed treatment comparisons. Statistics in Medicine 2004; 23:3105–3124. 5. Caldwell DM, Ades AE, Higgins JPT. Simultaneous comparison of multiple treatments: combining direct and indirect evidence. British Medical Journal 2005; 331:897–900. 6. Salanti G, Higgins JPT, Ades AE, Ioannidis JPA. Evaluation of networks of randomized trials. Statistical Methods in Medical Research 2008; 17(3):279–301. 7. Glenny AM, Altman DG, Song F, Sakarovitch C, Deeks JJ, D’Amico R, Bradburn M, Eastwood AJ, International Stroke Trial Collaborative Group. Indirect comparisons of competing interventions. Health Technology Assessment 2005; 9(26): 1–149. 8. Jansen JP, Fleurence R, Devine B, Itzler R, Barrett A, Hawkins N, Lee K, Boersma C, Annemans L, Cappelleri JC. Interpreting indirect treatment comparisons & network meta-analysis for health care decision-making: report of the ISPOR task force on good research practices - Part 1. Value in Health 2011; 14:417–428. 9. Hoaglin DC, Hawkins N, Jansen JP, Scott DA, Itzler R, Cappelleri JC, Boersma C, Thompson D, Larholt KM, Diaz M, Barret A. Conducting indirect treatment comparison and network meta-analysis studies: report of the ISPOR task force on indirect treatment comparisons good research practices - Part 2. Value in Health 2011; 14:429–437. 10. Dias S, Welton NJ, Sutton AJ, Ades AE. NICE DSU technical support document 2: a generalised linear modelling framework for pairwise and network meta-analysis of randomised controlled trials, 2011. Available from: http://www.nicedsu. org.uk [Accessed on 25 September 2012]. 11. Trikalinos TA, Olkin I. Meta-analysis of effect sizes reported at multiple time points: a multivariate approach. Clinical Trials 2012; 9(5):610–620. 12. Jones AP, Riley RD, Williamson PR, Whitehead A. Meta-analysis of individual patient data versus aggregate data from longitudinal clinical trials. Clinical Trials 2009; 6(1):16–27.

Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE 13. Ishak KJ, Platt RW, Joseph L, Hanley JA, Caro JJ. Meta-analysis of longitudinal studies. Clinical Trials 2007; 4(5): 525–539. 14. Peters JL, Mengersen KL. Meta-analysis of repeated measures study designs. Journal of Evaluation in Clinical Practice 2008; 14:941–950. 15. Wandel S, Jüni P, Tendal B, Nüesch E, Villiger PM, Welton NJ, Reichenbach S, Trelle S. Effects of glucosamine, chondroitin,or placebo in patients with osteoarthritis of hip or knee: network meta-analysis. British Medical Journal 2010; 341:c4675. 16. Dakin HA, Welton NJ, Ades AE, Collins S, Orme M, Kelly S. Mixed treatment comparison of repeated measurements of a continuous endpoint: an example using topical treatments for primary open-angle glaucoma and ocular hypertension. Statistics in Medicine 2011; 30(20):2511–2535. 17. Lu G, Ades AE, Sutton AJ, Cooper NJ, Briggs AH, Caldwell DM. Meta-analysis of mixed treatment comparisons at multiple follow-up times. Statistics in Medicine 2007; 26(20):3681–3699. 18. Royston P, Altman DG. Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling (with discussion). Applied Statistics 1994; 43:429–467. 19. Lambert PC, Smith LK, Jones DR, Botha JL. Additive and multiplicative covariate regression models for relative survival incorporating fractional polynomials for time-dependent effects. Statistics in Medicine 2005; 24:3871–3885. 20. Bossard N, Descotes F, Bremond AG, Bobin Y, De Saint Hilaire P, Golfier F, Awada A, Mathevet PM, Berrerd L, Barbier Y, Estève J. Keeping data continuous when analyzing the prognostic impact of a tumor marker: an example with cathepsin D in breast cancer. Breast Cancer Research and Treatment 2003; 82:47–59. 21. Berger U, Schafer J, Ulm K. Dynamic Cox modelling based on fractional polynomials: time-variations in gastric cancer prognosis. Statistics in Medicine 2003; 22:1163–1180. 22. Bagnardi V, Zambon A, Quatto P, Corrao G. Flexible meta-regression functions for modeling aggregate dose response data, with an application to alcohol and mortality. American Journal of Epidemiology 2004; 159:1077–1086. 23. Sauerbrei W, Royston P. Building multivariable prognostic and diagnostic models: transformation of the predictors by using fractional polynomials. Journal of the Royal Statistical Society Series A - Statistics in Society 1999; 162:71–94. 24. Sauerbrei W, Royston P, Look M. A new proposal for multivariable modelling of time-varying effects in survival data based on fractional polynomial time-transformation. Biometrical Journal 2007; 9:453–473. 25. Jansen JP. Network meta-analysis of survival data with fractional polynomials. BMC Medical Research Methodololog 2011; 11:61. 26. Jansen JP, Cope S. Meta-regression models to assess heterogeneity and inconsistency in network meta-analysis of survival outcomes. BMC Medical Research Methodololog 2012; 12:152. 27. Cooper NJ, Sutton AJ, Morris D, Ades AE, Welton NJ. Addressing between-study heterogeneity and inconsistency in mixed treatment comparisons: application to stroke prevention treatments in individuals with non-rheumatic atrial fibrillation. Statistics in Medicine 2009; 28:1861–1881. 28. Akaike H. Information theory and an extension of the maximum likelihood principle. Second International Symposium on Information Theory 1973; 1:267–281. 29. Dempster AP. The direct use of likelihood for significance testing. Statistics and Computing 1997; 7:247–252. 30. Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, Series B 2002; 64:583–639. 31. Dagenais S, Garbedian S. Systematic review of the prevalence of radiographic primary hip osteoarthritis. Clinical Orthopaedics and Related Research 2009; 467:623–637. 32. Pisoni C, Giardini A, Majani G, Maini M. International Classification of Functioning, Disability, and Health (ICF) core sets for osteoarthritis. A useful tool in the follow-up of patients after joint arthroplasty. European Journal of Physical and Rehabilitation Medicine 2008; 44:377–385. 33. Wenham CYJ, Conaghan PG. Imaging the painful osteoarthritic knee joint: what have we learned? Nature Clinical Practice Rheumatology 2009; 5:149–158. 34. Bellamy N, Campbell J, Robinson V, Gee T, Bourne R, Wells G. Viscosupplementation for the treatment of osteoarthritis of the knee. Cochrane Database of Systematic Reviews 2006; 19(2):CD005321. 35. Divine JG, Zazulak BT, Hewett TE. Viscosupplementation for knee osteoarthritis: a systematic review. Clinical Orthopaedics and Related Research 2007; 455:113–122. 36. Jordan KM, Arden NK, Doherty M, Bannwarth B, Bijlsma JW, Dieppe P, Gunther K, Hauselmann H, Herrero-Beaumont G, Kaklamanis P, Lohmander S, Leeb B, Lequesne M, Mazieres B, Martin-Mola E, Pavelka K, Pendleton A, Punzi L, Serni U, Swoboda B, Verbruggen G, Zimmerman-Gorska I, Dougados M, Standing Committee for International Clinical Studies Including Therapeutic Trials ESCISIT. EULAR recommendations 2003: an evidence based approach to the management of knee osteoarthritis: report of a task force of the Standing Committee for International Clinical Studies Including Therapeutic Trials (ESCISIT). Annals of the Rheumatic Diseases 2003; 62:1145–1155. 37. Simon LS, Lipman AG, Jacox AK, Caudill-Slosberg M, Gill LH, Keefe FJ, Kerr KL, Minor MA, Sherry DD, Vallerand AH, Vasudevan S. Guidelines for the management of pain in osteoarthritis, rheumatoid arthritis, and juvenile chronic arthritis. In APS Clinical Practice Guidelines Series, (2nd edn). American Pain Society: Glenview, IL; 2002. 38. Zhang W, Moskowitz RW, Nuki G, Abramson S, Altman RD, Arden N, Bierma-Zeinstra S, Brandt KD, Croft P, Doherty M, Dougados M, Hochberg M, Hunter DJ, Kwoh K, Lohmander LS, Tugwell P. OARSI recommendations for the management of hip and knee osteoarthritis, part II: OARSI evidence-based, expert consensus guidelines. Osteoarthritis Cartilage 2008; 16:137–162. 39. Altman RD, Moskowitz R. Intraarticular sodium hyaluronate (Hyalgan) in the treatment of patients with osteoarthritis of the knee: a randomized clinical trial. Hyalgan study group. The Journal of Rheumatology 1998; 25(11):2203–2212. 40. Altman RD, Rosen JE, Bloch DA, Hatoum HT, Korner P. A double-blind, randomized, saline-controlled study of the efficacy and safety of EUFLEXXA for treatment of painful osteoarthritis of the knee, with an open-label safety extension (the FLEXX trial). Seminars in Arthritis and Rheumatism 2009; 39(1):1–9.

Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015

J. P. JANSEN, M. C. VIEIRA AND S. COPE 41. Carrabba M, Paresce E, Angelini M, Re KA, Torchiana EEM, Perbellini A. The safety and efficacy of different dose schedules of hyaluronic acid in the treatment of painful osteoarthritis of the knee with joint effusion. European Journal of Rheumatology & Inflammation 1995; 15(1):25–31. 42. Corrado EM, Peluso GF, Gigliotti S, De Durante C, Palmieri D, Savoia N, Oriani GO, Tajana GF. The effects of intra-articular administration of hyaluronic acid on osteoarthritis of the knee: a clinical study with immunological and biochemical evaluations. European Journal of Rheumatology & Inflammation 1995; 15(1):47–56. 43. Cubukçu D, Ardiç F, Karabulut N, Topuz O. Hylan G-F 20 efficacy on articular cartilage quality in patients with knee osteoarthritis: clinical and MRI assessment. Clinical Rheumatology 2005; 24(4):336–341. 44. Dickson DJ, Hosie G, English JR. A double-blind, placebo controlled comparison of hylan G-F 20 against diclofenac in knee osteoarthritis. Journal of Clinical Research 2001; 4:41–52. 45. Dougados M, Nguyen M, Listrat V, Amor B. High molecular weight sodium hyaluronate (hyalectin) in osteoarthritis of the knee: a 1 year placebo-controlled trial. Osteoarthritis Cartilage 1993; 1(2):97–103. 46. Henderson EB, Smith EC, Pegley F, Blake DR. Intra-articular injections of 750 kD hyaluronan in the treatment of osteoarthritis: a randomised single centre double-blind placebo-controlled trial of 91 patients demonstrating lack of efficacy. Annals of the Rheumatic Diseases 1994; 53(8):529–534. 47. Huskisson EC, Donnelly S. Hyaluronic acid in the treatment of osteoarthritis of the knee. Rheumatology (Oxford) 1999; 38(7):602–607. 48. Jubb RW, Piva S, Beinat L, Dacre J, Gishen P. A one-year, randomised, placebo (saline) controlled clinical trial of 500-730 kDa sodium hyaluronate (Hyalgan) on the radiological change in osteoarthritis of the knee. International Journal of Clinical Practice 2003; 57(6):467–474. 49. Karlsson J, Sjogren LS, Lohmander LS. Comparison of two hyaluronan drugs and placebo in patients with knee osteoarthritis. A controlled, randomized, double-blind, parallel-design multicentre study. Rheumatology (Oxford) 2002; 41(11): 1240–1128. 50. Kirchner M, Marshall D. A double-blind randomized controlled trial comparing alternate forms of high molecular weight hyaluronan for the treatment of osteoarthritis of the knee. Osteoarthritis and Cartilage 2006; 14:154–162. 51. Scale D, Wobig M, Wolpert W. Viscosupplementation of osteoarthritic knees with hylan: a treatment schedule study. Current Therapeutic Research, Clinical and Experimental 1994; 55(3):220–232. 52. Stitik TP, Blacksin MF, Stiskal DM, Kim JH, Foye PM, Schoenherr L, Choi ES, Chen B, Saunders HJ, Nadler SF. Efficacy and safety of hyaluronan treatment in combination therapy with home exercise for knee osteoarthritis pain. Archives of Physical Medicine and Rehabilitation 2007; 88(2):135–141. 53. Wobig M, Dickhut A, Maier R, Vetter G. Viscosupplementation with hylan G-F 20: a 26-week controlled trial of efficacy and safety in the osteoarthritic knee. Clinical Therapeutics 1998; 20(3):410–423. 54. Wobig M, Bach G, Beks P, Dickhut A, Runzheimer J, Schwieger G, Vetter G, Balazs E. The role of elastoviscosity in the efficacy of viscosupplementation for osteoarthritis of the knee: a comparison of hylan G-F 20 and a lower-molecular-weight hyaluronan. Clinical Therapeutics 1999; 21(9):1549–1562. 55. Lunn D, Spiegelhalter D, Thomas A, Best N. The BUGS project: evolution, critique, and future directions. Statistics in Medicine 2009; 28:3049–3067. 56. Dias S, Welton NJ, Sutton AJ, Caldwell DM, Guobing L, Ades AE. NICE DSU technical support document 4: inconsistency in networks of evidence based on randomised controlled trials, 2011. Available from: http://www.nicedsu.org.uk [Accessed on 25 September 2012]. 57. Riley RD. Multivariate meta-analysis: the effect of ignoring within-study correlation. Journal of the Royal Statistical Society Series A - Statistics in Society 2009; 172:789–811.

Copyright © 2015 John Wiley & Sons, Ltd.

Statist. Med. 2015