Bayesian one‐step IPD network meta‐analysis of time‐to‐event data using Royston‐Parmar models

Network meta‐analysis (NMA) combines direct and indirect evidence from trials to calculate and rank treatment estimates. While modelling approaches for continuous and binary outcomes are relatively well developed, less work has been done with time‐to‐event outcomes. Such outcomes are usually analysed using Cox proportional hazard (PH) models. However, in oncology with longer follow‐up time, and time‐dependent effects of targeted treatments, this may no longer be appropriate. Network meta‐analysis conducted in the Bayesian setting has been increasing in popularity. However, fitting the Cox model is computationally intensive, making it unsuitable for many datasets. Royston‐Parmar models are a flexible alternative that can accommodate time‐dependent effects. Motivated by individual participant data (IPD) from 37 cervical cancer trials (5922 women) comparing surgery, radiotherapy, and chemotherapy, this paper develops an IPD Royston‐Parmar Bayesian NMA model for overall survival. We give WinBUGS code for the model. We show how including a treatment‐ln(time) interaction can be used to conduct a global test for PH, illustrate how to test for consistency of direct and indirect evidence, and assess within‐design heterogeneity. Our approach provides a computationally practical, flexible Bayesian approach to NMA of IPD survival data, which readily extends to include additional complexities, such as non‐PH, increasingly found in oncology trials.

A. Worked example using the cervical cancer network Treatment parameterisation For the cervical cancer network, let trt1 i be the indicator variable for RT v CTRT, trt2 i be the indicator variable for RT v CT+RT and trt3 i be the indicator variable for RT v CT+S. The treatment contrast variables are defined as: 1 if patient was randomised to CTRT and is from a trial comparing RT and CTRT 0 otherwise if patient was randomised to CT+RT and is from a trial comparing RT  Stata code for initial data processing In this subsection we present Stata code for the initial data processing steps which includes loading data, setting data as survival data and generating the treatment contrast variables.

Calculating basis functions
In this subsection we explain how to calculate the basis functions introduced in Section 3.1.
For a patient i with survival (or censoring time) ln(t i ) and a spline model with p interior knots where the location of the knots are k 0 , k 1 , . . . , k p , k p+1 (k 0 < k 1 < ... < k p < k p+1 ) and k 0 and k p+1 are the boundary knots and where (x) + equals x if x > 0 and 0 otherwise, the basis functions v 0 ln(t i ) , v 1 ln(t i ) , v 2 ln(t i ) , . . . , v p−1 ln(t i ) , v p ln(t i ) are calculated as follows: into u 0 ln(t i ) , u 1 ln(t i ) , . . . , u p−1 ln(t i ) , u p ln(t i ) which can then be used in equations (1), (2),

Stata code for calculating the basis functions
This subsection includes the Stata code for calculating the basis functions. We start by first specifying the location of the internal knots. In this example code we assume that all trials have the same two internal knot locations (based on percentiles of survival time). If the knot locations are chosen individually for each trial (as in the cervical cancer network) then these should be specified individually for each trial.
*** Set up location of knots for each trial } * Here we assume all trials have the same two internal knot locations * at the 33rd and 67th percentiles of uncensored survival times (0 33 67 100) replace

Stata code for running WinBUGS
In the next section of Stata code we use the winbugs suite of commands to prepare the data for Win-BUGS, write a script file (which will tell WinBUGS which files to run), run the model in WinBUGS and load the results back in Stata.

WinBUGS model code
This next section specifies the one-step RTE NMA model (6) that was run in WinBUGS using the cervical cancer network. In the model i is patient and j is trial. This model specification also shows how the spline can be extended to fit models (8), (9) and (10). Note for these models the data file will need to be modified to include the treatment-ln(time) interactions or the indicator variable for the inconsistency parameter.  (10) # Consider the treatment loop in Figure 1 to be coded as: # Extension to (7) #  LGOG -0.44 0% 0.6 68% n/a n/a 1.16 100% 15