Abstract
Multiple randomized controlled trials, each comparing a subset of competing interventions, can be synthesized by means of a network meta-analysis to estimate relative treatment effects between all interventions in the evidence base. Often there is an interest in estimating the relative treatment effects regarding time-to-event outcomes. Cancer treatment effectiveness is frequently quantified by analyzing overall survival (OS) and progression-free survival (PFS). In this paper we introduce a method for the joint network meta-analysis of PFS and OS that is based on a time-inhomogeneous tri-state (stable, progression, and death) Markov model where time-varying transition rates and relative treatment effects are modeled with known parametric survival functions or fractional polynomials. The data needed to run these analyses can be extracted directly from published survival curves. We demonstrate use by applying the methodology to a network of trials for the treatment of non-small-cell lung cancer. The proposed approach allows the joint synthesis of OS and PFS, relaxes the proportional hazards assumption, extends to a network of more than two treatments, and simplifies the parameterization of decision and cost-effectiveness analyses.
1 Introduction
Randomized controlled trials (RCTs) are considered the most appropriate study design to obtain evidence regarding the relative treatment effects of an intervention. However, an individual RCT rarely includes all alternative interventions of interest, and as such does not provide all the information needed to select the best option among alternatives. Typically, the evidence base consists of multiple RCTs where each of the available studies compares a subset of the alternative interventions of interest. If each of these trials has at least one intervention in common with another trial such that the evidence base is represented by a single connected network, a network meta-analysis (NMA) can provide relative treatment effects between all the competing interventions in the evidence base ^{1}.
Often there is an interest in estimating the relative treatment effects of alternative interventions regarding time-to-event outcomes. For example in oncology, treatment efficacy is often quantified by analyzing time from treatment initiation to the occurrence of a particular event. Very commonly studies report data on overall survival (OS), where the event is death from any cause, and on progression-free survival (PFS), where the event is death from any cause or disease progression, whichever occurred first.
NMA of time-to-event outcomes with a single effect measure per study are based on the proportion of patients alive at a specific time point, median survival, or reported hazard ratio (HR) ^{2}. The limitation of a NMA of survival at a specific time point is that we only focus on the cumulative effect of treatment at that time point and ignore the variation in effects over time up to, as well as beyond, that time point. NMAs of median survival times have similar limitations. The HR summarizes the treatment effect for the complete follow-up period of the trials, but only represents the treatment effect for each time point if the proportional hazards (PH) assumption holds. If the PH assumption is violated, trial specific HRs represent an average effect over the follow-up period, which can cause biased estimates in a NMA if trials have different lengths of follow-up.
As an alternative to a NMA with a univariate treatment effect measure, we can also use a multivariate treatment effect measure that describes how the relative treatment effects change over time ^{2}. Ouwens et al., Jansen, and Cope et al presented methods for NMA of time-to-event outcomes where the hazard functions of the interventions in a trial are modeled using known parametric survival functions or fractional polynomials and the difference in the parameters are considered the multi-dimensional treatment effects, which are synthesized across studies ^{3–6}. By incorporating time-related parameters, these NMA models can be fitted more closely to the available data.
Both PFS and OS of an intervention determine its value and can inform decision-making. In combination with a baseline survival function for a reference treatment, the multivariate NMA models embedded in parametric survival functions can form the basis for partitioned-survival cost-effectiveness models. Frequently, the pooled PFS and OS curves need to be extrapolated over time in order to obtain estimates of the expected quality adjusted life-years before and after disease progression. Since the separate meta-analyses of PFS and of OS data ignore the correlation between the outcomes, any required extrapolation may result in possible crossing of PFS and OS curves. A state-transition model with three health states – stable (pre-progression), progression, and death – with parametric hazard functions for the three corresponding transisitions avoids this issue. If we have individual patient data (IPD) regarding time to progression, time to death, and censoring for all trials included in the NMA then we can estimate these hazard functions using a statistical model with the same tri-state structure and avoid any inconsistency between the clinical evidence synthesis and the economic evaluation ^{7}. Reality though is that for most, if not all, trials there is no access to IPD and the synthesis has to be based on reported summary findings. Although reported Kaplan-Meier curves for PFS and OS can be digitized and a dataset of “virtual” IPD event-times can be created with the algorithm by Guyot et al., it does not provide the information needed to determine which time-to-progression data point corresponds to which to time-to-death data point^{8}. Markov-state-transition NMA models have been presented for disease progression^{9,10} and competing risks ^{11} based on aggregate level data, but these models assumed constant hazards for transitions between states.
In this paper we introduce a method for the joint NMA of PFS and OS that is based on a tri-state (stable, progression, and death) transition model where time-varying hazard rates and relative treatment effects are modeled with known parametric survival functions or fractional polynomials. We illustrate parameter estimation based on aggregate level data.
2 Multi-state network meta-analysis framework
At any time u, patients in study i randomized to treatment arm k can be in one of three health states: alive with stable disease (i.e. not progressed), alive and progressed, and dead expressed with S_{ik}(u), P_{ik}(u) and D_{ik}(u) respectively and illustrated in Figure 1 , and are the hazard rates for disease progression, dying post-progression, and dying pre-progression.
A multi-state NMA that explicitly estimates each possible transition in a tri-state model and modeling time-varying hazard rates and relative treatment effects with known parametric survival functions or fractional polynomials can be expressed as follows: where p_{1}, …, p_{B} are fractional powers and the round bracket notation denotes the Box-Tidwell transformation: u^{(p)} = u^{p} if p ≠ 0 and u^{(p)} = ln(u) if p = 0. Equation 1 also includes the situation of repeated powers, where p_{x} = p_{y} for at least 1 pair of indices (x, y), 1 ≤ x < y ≤ B. In this situation, is used instead of itself. A complete set of flexible fractional polynomials can be created with {(p_{1}, …, p_{B−1}) = (−2, −1, −0.5, 0, 0.5, 1, 2).
α_{1,ik}, α_{2,ik}, …, α_{a,ik} are regression coefficients that represent the scale and shape parameters of the log hazard function describing the stable-to-progression transition in study i for treatment arm k. α_{a+1,ik}, α_{a+2,ik}, …, α_{b,ik} are the regression coefficients that represent the log-hazard function for the stable-to-death transition. α_{b+1,ik}, α_{b+2,ik}, …, α_{B,ik} are the regression coefficients that represent the scale and shape parameters of the log hazard function describing the progression-to-death in study i for treatment arm k.
The µ_{·i} reflect the study effects regarding the scale and shape parameters in each study i. The δ_{·,ik} are the study specific true underlying relative treatment effects for the treatment in arm k relative to the treatment in arm 1 of that trial (with δ_{·,i1} = 0) regarding the scale and shape of the log hazard function for the different transitions, which are drawn from a normal distribution with the mean effect for treatment t expressed in terms of the overall reference treatment 1,, and with a between-study-heterogeneity covariance matrix Σ. A fixed effects model is obtained by replacing with .
3 Illustrative example
3.1 Evidence base
An example of the multi-state models is presented for a NMA of first line treatment of adult patients with metastatic EGFR+ non-squamous non-small-cell lung cancer (NSCLC) with gefitinib, erlotinib, afatinib, dacomitinib, or platinum-based doublet chemotherapy regimens. Thirteen RCTs were obtained with a systematic literature review (ARCHER1050^{12,13}; LUX-LUNG 7^{14,15}; LUX-LUNG 3^{16,17}; LUX-LUNG 6^{17,18}; EURTAC ^{19–21}; ENSURE^{22}; OPTIMAL^{23,24}; First-SIGNAL ^{25}; WJTOG3405^{26}; IPASS ^{27,28}; NEJ002^{29,30}; Han2017^{31}; Yang2014^{32,33}). The evidence network is presented in Figure 2 and the trial-specific PFS and OS curves are provided as supplementary information in Section C.1.
A NMA provides estimates of relative treatment effects between the competing interventions (i.e. hazard ratios). In order to obtain estimates for the transition rates over time between health states for each treatment, we first need to obtain these estimates for an overall reference treatment, defined as treatment 1. Next, the hazard ratios of each treatment relative to treatment 1 obtained with the NMA are applied to these transition rates for treatment 1. As a final step, these time-varying transition rates can be transformed into the distribution S, P, and D over time, and PFS and OS curves. In this example gefitinib is defined as treatment 1, and a meta-analysis of gefitinib arms was performed as well. In the following sections we present the NMA model as well as the meta-analysis model used for this example.
3.2 Network meta-analysis
3.2.1 Model
The following model was used for the NMAs, which is a simplification of Equation 1 to faciliate parameter estimation: where u^{0} = ln(u) and d_{1,11} = 0, d_{2,11} = 0, d_{3,11} = 0, and d_{4,11} = 0
When p_{1} = 0 and α_{3,ik} = 0, the log-hazard functions follow a Weibull distribution. When p_{1} = 1 and α_{3,ik} = 0 these log-hazard functions follow a Gompertz distribution. When {(p_{1}, p_{2}) = (0, 0), (0, 1)} and α_{3,ik} ≠ 0, the log-hazard functions follow a second order polynomial that are extensions of the Weibull and Gompertz model to allow for arc- and bathtub shaped log-hazard functions. When α_{6,ik} = 0 this transition follows an exponential distribution. When α_{6,ik} ≠ 0 the transition rates over time are modeled with a first order fractional polynomial of which Weibull and Gompertz are special cases.
With this model we assume there is one between-study heterogeneity parameter related to the relative treatment effects for . The δ_{1,ik} are drawn from a normal distribution with the mean effect for treatment t expressed in terms of the overall reference treatment 1,, and between study heterogeneity . The treatment specific relative effects regarding the first shape parameter of the log hazard function for the stable-to-progression transition was assumed to be fixed. Incorporating the treatment specific relative effects regarding the second shape parameter frequently results in unstable parameter estimates and was removed from the model. The same can be argued for the treatment effect regarding the shape parameter of the log hazard function for the progression-to-death transition. If it is assumed that treatment has only a direct effect on the transitions from stable to progression, which can reasonably be defended when a particular treatment upon disease progression is discontinued, the model can be further simplified by setting .
The random effects model presented with Equation 1 does not account for correlation between trial-specific δ_{1,ik}s in multiple-arm trials (>2 treatments). A random effects model with only a heterogeneity parameter for can be easily extended to fit trials with three or more treatment arms by decomposition of a multivariate normal distribution as a series of conditional univariate distributions^{34}: The vector of random effects follows a multivariate normal distribution, a_{i} represents the number of arms in trial i (a_{i} = 2, 3,) and . The conditional univariate distribution for the random effect of arm k > 2, given all arms from 2 to k â1, is
3.2.2 Likelihood
The model parameters were estimated based on the conditional survival probabilities regarding PFS and OS obtained from the published KM curves. The total follow-up time can be partitioned into M successive non-overlapping intervals indexed by m = 1, …, M. We refer to interval m as U_{m} and write u ∈ U_{m} to denote u_{m} ≤ u < u_{m+1}. The length of U_{m} is Δu_{m} = u_{m+1} − u_{m}. For each interval m, a binomial likelihood for the conditional survival probabilities regarding PFS and OS at time point u relative to the time point at the beginning of the interval u_{m} can be described according to: where are the observed number of patients who have not yet experienced progression or death at time u in the m^{th} interval in study i for treatment arm k and are the observed number of patients who have not died at time u in that interval. is the underlying conditional survival probability regarding PFS and is the underlying conditional survival probability regarding OS, and are the corresponding sample sizes at the beginning of the interval. A description how to create interval data based on KM curves is provided in Appendix A).
For the m^{th} interval, the conditional probabilities and are related to the proportion of patients who are progression free (stable disease) S_{ik}(u) and the proportion of patients with progressed disease P_{ik}(u) according to: Arbitrary hazard functions can be approximated with a set of discontinuous constant hazard rates over relative short successive time intervals. For each interval m, S_{ik}(u), P_{ik}(u), and death D_{ik}(u) are related to the hazards according to the following set of differential equations (See Appendix B): In order to estimate the three parameters , and for each interval m, we need to define Equation 5, Equation 6 and Equation 7 for at least two time points per interval. We use the mid point , which we define as , and the time point at the end of the interval u_{m+1}. The obtained estimates of the hazards for interval m were assigned to the time point for Equation 2.
3.2.3 Prior distributions
The following prior distributions for the parameters of the model expressed with Equation 1 were used:
3.3 Meta-analysis of absolute effects with overall reference treatment
The following fixed effects model was used to estimate transition rates for gefitinib: We used the same data structure, likelihood, and link functions as used for the NMA. (See Equation 5,Equation 6, and Equation 7). The prior distribution for this model was:
3.4 Parameter estimation
The parameters of the different models were estimated using a Markov Chain Monte Carlo (MCMC) method implemented in the JAGS software package^{35}. All JAGS analyses were run using R statistical software ^{36}. See Section C.2 for the JAGS code for one of the models used to estimate relative treatment effects.
The residual deviance and the deviance information criterion (DIC) were used to compare the goodness-of-fit of the competing models. The DIC provides a measure of model fit that penalizes model complexity. In general, a more complex model results in a better fit to the data, demonstrating a smaller residual deviance. The model with the better trade-off between fit and parsimony has a lower DIC. A difference in the DIC of about 10 points can be considered meaningful.
3.5 Results
For this example we used the following competing models for the analysis of treatment 1 (gefitinib):
- SP 2nd order FP(01); SD exponential; PD Weibull. This is a fixed effects 2nd order fractional polynomial with time transformations according to p_{1} = 0 and p_{2} = 1 for the stable-to-progression (SP) transition; an exponential distribution for the stable-to-death (SD) transition; and a Weibull distribution for the progression-to-death (PD) transition. (See Equation 9)
- SP 2nd order FP(01); SD exponential; PD exponential
- SP 2nd order FP(00); SD exponential; PD Weibull (See Equation 9 with p_{1} = p_{2} = 0)
- SP Weibull; SD exponential; PD Weibull
- SP Weibull; SD exponential; PD exponential
- SP Gompertz; SD exponential; PD Weibull
The model fit statistics are presented in Table 1. With the exception of the model where a Gompertz distribution was used, the models for the meta-analysis have a similar deviance and DIC. The parameter estimates of four competing models that show the greatest variation in time-varying hazards between the three health states are presented in Table 2. The actual time-varying hazard rates and corresponding PFS and OS curves are plotted in Figure 3.
For the NMA the following models were evaluated:
- SP 2nd order FP(01) FE3; SD exponential; PD Weibull FE1(scale). This is a 2nd order fractional polynomial with time transformations according to p_{1} = 0 and p_{2} = 1 for the stable-to-progression transition with fixed effects relative treatment effects on all 3 parameters (α_{1}, α_{2}, α_{3}); an exponential distribution for the stable-to-death transition without relative treatment effects; and a Weibull distribution for the progression-to-death transition with fixed effects relative treatment effects only on the scale parameter (α_{5}).
- SP 2nd order FP(01) FE3; SD exponential; PD Weibull. This model only has relative treatment effects applied to the stable-to-progression transition
- SP 2nd order FP(00) FE3; SD exponential; PD exponential
- SP Weibull FE2; SD exponential; PD Weibull FE1(scale)
- SP Weibull FE2; SD exponential; PD Weibull
- SP Weibull FE2; SD exponential; PD exponential
- SP 2nd order FP(01) RE3; SD exponential; PD Weibull FE1(scale). This is the random effects model according to Equation 2 with p_{1} = 0 and p_{2} = 1.
- SP Weibull RE2; SD exponential; PD Weibull FE1(scale)
The NMA models that assumed a 2nd order fractional polynomial for the stable-to-progression transitions had a lower deviance than the simpler models assuming a Weibull function, but cannot be considered a meaninful improvement when factoring in model complexity, as indicated by the DIC. Comparing models that assumed a relative treatment effect for the progression-to-death transition with the corresponding models without indicates that a relative effect is an important component to include for this transition. The 2nd order fractional polynomial and Weibull random effects models performed better than their fixed effects equivalents, indicating that incorporating between-study heterogeneity is important.
Parameter estimates for the random effects 2nd order fractional polynomial model (SP 2nd order FP(01) RE3; SD exponential; PD Weibull FE1(scale)), the random effects Weibull model (SP Weibull RE2; SD exponential; PD Weibull FE1(scale), the corresponding fixed effects model (SP Weibull FE2; SD exponential; PD Weibull FE1(scale)), and the simplest model (SP Weibull FE2; SD exponential; PD exponential) are presented in Table 3. The corresponding time-varying HRs with treatment relative to treatment 1 for the stable to progression transitions are presented in Figure 4 and the constant HRs for the progression-to-death transition in Figure 5. (Please note that we did not assume a relative treatment effect for the stable-to-death transition). Applying the relative treatment effect parameters describing the HRs over time obtained with these NMA models to the corresponding models for the analysis of treatment 1 we obtain the PFS and OS curves by treatment, as presented in Figure 6. In order to illustrate the width of the 95 credible intervals of these survival curves due to the uncertainty in the time-varying HRs, we ignored the uncertainty for the reference treatment 1. (If the uncertainty of the meta-analysis would have been incorporated as well, the 95 credible intervals would have been a bit wider.)
4 Discussion
With this paper we present a method for the joint NMA of PFS and OS that is based on a tri-state transition model. This method extends existing parametric NMA methods for time-to-event data ^{3–6} by defining the structural relationship between PFS and OS according to the stable, progression, and death states that define the course of disease over time. Instead of modeling the time-varying hazard rates for PFS and OS separately, we model the time-varying transition rates between the three health states simultaneously. The primary advantage of this evidence synthesis framework is that estimates for PFS and OS remain consistent over time.
For the illustrative example we used a subset of the theoretically possible competing models available within the proposed framework by making simplifying structural assumptions to facilitate parameter estimation. First, we used an exponential distribution for the stable-to-death transition, informed by the notion that this rate is representative of general mortality which is low relative to disease progression related mortality and constant over time given the life-expectancy of this population. Second, we assumed that a relative treatment effect for this transition was not needed reflecting the belief that differences in survival between treatments are only due to differences in delayed or avoided tumor progression, and differences in adverse event rates do not influence mortality. Third, we assumed that a function more complex than a first order fractional polynomial with a constant relative treatment effect was not required for the progression-to-death transition given that this transition is conditional upon experiencing progression and modeled in relation to follow-up time. One could argue that a relative treatment effect for this transition is not needed when treatment is discontinued upon progression. However, the DIC indicated that adding this parameter to the models resulted in a meaningful improvement for the example analyses, which is primarily related to the PBDC trials. A reason to include a treatment effect parameter for the progression-to-death transition is treatment cross-over upon progression in a subset of trials or if post-progression treatment between trial arms differ. Fourth, out of the possible fractional polynomials, we only evaluated exponential, Weibull and Gompertz models and their extensions where the additional parameter related the log hazard to time or log-time. We did not consider any of the negative power transformations of time, primarily because these functions do not link to known survival distributions and the 2nd order models we did use have already the flexibility to capture arc-shaped hazard functions.
This brings us to the point of model selection. Evaluating the fit of all possible competing models available within the proposed framework to the data is not feasible. Future research is needed to inform a model selection strategy or algorithm that results in a set of models that is likely to cover the distribution of transition rates between the health states, results in realistic extrapolations over time, and, given the computational burden of the more complex models for large datasets, can be evaluated in a reasonable amount of time.
The primary reason to propose the method described in this paper is to facilitate parameterization of multi-state cost-effectiveness models. More specifically, the proposed evidence synthesis framework relates directly to clock forward time-inhomogeneous Markov models where treatment specific transition rates between health states are only a function of time in the model. A frequently used approach for cost-effectiveness analysis of cancer treatments are partitioned survival models. However, the main limitation is that extrapolated parametric PFS and OS curves for a given treatment may cross. This will not be case with Markov models and, as such, are preferred as long as time-varying transition rates between health states can be estimated that reflect the actual PFS and OS of the treatments compared. As far as we know, the method presented in this paper is the first to facilitate this based on reported aggregate level data. In order to obtain the parameter estimates for a cost-effectiveness analysis we need to define a baseline model and a NMA model. The baseline model provides estimates for the absolute effect with the reference treatment, which in this case are the time-varying log-hazard rates between each of the three health states. The NMA model provides estimates of the relative treatment effects of each intervention in the network relative to the reference treatment, which in this case are the time-varying log hazards ratios. The absolute effect with each treatment is obtained by adding the relative treatment effects from the NMA to the absolute effect with the reference treatment from the baseline model, and subsequently transforming these to the natural scale by inverting the log-link function ^{37}. In the current example we used the RCT evidence base to estimate the baseline meta-analysis model as well as the NMA model. However, for an actual cost-effectiveness analysis it may be more appropriate to use a different evidence base that better reflects expected outcomes with the reference treatment for the baseline model, such as a long-term routine practice observational studies.
The estimates obtained with the proposed evidence synthesis models can also be used in semi-Markov individual simulation models (i.e. models where some transitions are affected by time in an intermediate state). However, for each of the transitions for which the “clock is reset” a separate multi-state (network) meta-analysis needs to be performed. For example, imagine a cost-effectiveness model of sequential cancer treatment consisting of four health states: 1) stable disease with first line treatment, 2) progression with first line treatment/stable disease with second line treatment, 3) progression with second line treatment, and 4) death. First-line treatment transitions from stable-to-progression and stable-to-death are estimated with one multi-state NMA model based on first line trials, and the second-line treatment transitions from stable-to-progression and stable-to-death are estimated with another multi-state NMA model based on second line trials.
All studies provided PFS and OS Kaplan-Meier data in the example analyses. In principle, the NMA model can be extended to create a shared-parameter model to incorporate studies that only provide information for PFS or only for OS. Studies with only PFS data provide evidence regarding the stable-to-progression and stable-to-death transitions and contribute to estimating the corresponding treatment specific hazard ratios if these are assumed fixed or exchangeable across all studies providing direct or indirect evidence for that particular intervention. (When a meta-analysis of absolute effects with the overall reference treatment is performed, the fixed effects or exchangeability assumption applies to the transition rates.) Incorporating studies that only provide OS data for a particular intervention in the NMA will require the additional assumption of fixed or exchangeable rates for one of the transitions across all studies for that intervention, if treatment is assumed to impact more than just the stable-to-progression transition in order to facilitate parameter estimation. A related topic for future research is whether and how this framework can be used to validate PFS as a surrogate for OS and to predict OS for novel interventions for which only mature PFS is available. This will be of great benefit for cost-effectiveness analyses.
5 Conclusion
We introduced a method for the joint meta-analysis of PFS and OS that is based on a non-homogenous Markovian tri-state transition model. Arbitrary hazard rate functions can be approximated by piecewise constant hazard rates at successive time intervals, and are flexibly modeled as (fractional) polynomial functions of time. The proposed approach relaxes the proportional hazards assumption, extends to a network of more than two treatments, and simplifies the parameterization of decision and cost-effectiveness analyses. The data needed to run these analyses can be extracted directly from published survival curves.
Data Availability
Dataset is available upon reasonable request.
6 Conflict of interest statement
All authors have completed the ICMJE uniform disclosure form at www.icmje.org/coi_disclosure.pdf and declare: no financial support from any organisation for the submitted work; JPJ is a part-time employee and shareholder of the Precision Medicine Group (Precision HEOR); DI is currently a full time employee of Genentech, the research was primarily performed while he was an employee of the Precision Medicine Group (Precision HEOR); TT has nothing to disclose.
7 Funding
The authors did not receive any funding related to the work presented.
Appendices
A Constructing dataset for analyses
Data inputs required are the coordinates extracted from digitally scanned PFS and OS Kaplan-Meier curves: time points (u),corresponding survival probabilities (S(u)), and corresponding population size at risk (n_{u}). These points must capture all steps in the curve, and may require adjustments to the extracted coordinates to ensure the survival probabilities are decreasing with time. For both curves it should include the times at which numbers at risk are reported below the curve.
The total follow-up time can be partitioned into M successive non-overlapping intervals indexed by m = 1, …, M. We refer to interval m as U_{m} and write u ∈ U_{m} to denote u_{m} ≤ u < u_{m+1}. The length of U_{m} is Δu_{m} = u_{m+1} − u_{m}. For each time interval m, we want to obtain three data points: At the beginning of the interval, u_{m}; at the mid point , which we define as ; and at the end of the interval, u_{m+1}. It is desirable to have the time intervals defined in such a way that (some) of these time points are aligned with the time point for which the size of the at-risk population is reported below the published Kaplan-Meier curves, and are the same for PFS and OS where available. For the current study, we used intervals with a length of 2 months.
If no PFS or OS proportion have been recorded for a specific time point of interest (i.e. whole months), a corresponding value for S(u) can be obtained by linear interpolation of the first available extracted scanned survival proportions before and after this time point.
When the population n_{u} is not reported below the PFS and OS Kaplan-Meier curve for certain time points u, it can be imputed. First, based on the reported size of the at-risk population at subsequent time points (n_{u+1}), n_{u} will be estimated according to . With this ‘backward calculation’ approach we implicitly assume that censoring occurs before the events happen within a time interval. However, this approach is not feasible if there is no information regarding the at-risk population for time intervals beyond the at-risk population reported at a certain time point. In other words, this approach is only feasibly for intervals up to the latest time point for which population is reported. Next, n_{u} will be estimated according to . The disadvantage of this ‘forward calculation’ approach is that censoring is ignored and the sample size potentially too large for those timepoints. For intervals where both and was calculated, the actual estimate for the population at-risk is calculated as: to ensure the sample size is not overestimated. For time points where could not be calculated,.
Based on the subsequent S(u) for the three points at each interval (i.e. S(u_{m}),, and S(u_{m+1})), two conditional survival proportions are obtained: and . The corresponding sample sizes are defined as . The corresponding observed number of patients who have not yet experienced progression or death are calculated according to .
Applying this algorithm to PFS and OS of each arm i of each trial k, we get a data set with and .
We set-up the event dataset such that every row represents one time interval with and corresponding to and u_{m+1}. In addition, each row has a variable related to follow-up time , two variables related to and u_{m+1} − u _{m} the study number, and study-arm number within that study.
In additon to the event dataset, we create a study dataset indicating the compared interventions in each study along with the number of study arms.
B States and between-state transition rates
B.1 Dynamic transitions - Problem specification
Figure 1 represents a closed dynamic system (S_{ik}(u) + P_{ik}(u) + D_{ik}(u) = 1) whose evolution is determined by a known initial condition at time u = 0 and three differential equations: With , and the time-varying hazard rates for the transitions in the figure.
B.2 Aproximating arbitrary , and
We can approximate arbitrary hazard rate functions with a set of discontinuous constant hazard rates over successive time intervals. We prefer this approximation because the system Equation A1 can be solved analytically when the transition rates are constant using the the eigenvalue method for first-order differential equations. For u ∈ U_{m} Equation A1 become:
B.3 Analytic solutions for S_{ik}(u), P_{ik}(u), and D_{ik}(u) where u ∈ U_{m}
Write the system in Equation A2 in matrix form: with the obvious notational correspondence between the two equations. For u ∈ U_{m} the system is homogenous and its general solution is the superposition: where λ_{1,ik}, λ_{2,ik}, and λ_{3,ik} are the eigenvalues of the coefficient matrix A_{ik}.v_{1ik}, v_{2ik} and v_{3ik} are the corresponding eigenvectors, and c_{1,ik}, c_{2,ik}, c_{3,ik} scalar constants to be identified from the initial condition in Equation A2. In our case: The eigenvectors are:
B.3.1 Identification of constants in the general solution
The constants c_{1,ik}, c_{2,ik}, c_{3,ik} are identied from the proportions at the beginning of U_{m}. Setting u = u_{m} in the general solution, and using the initial condition in Equation A1 and Equation A2 we obtain: where v_{xy,ik} is element x of eigenvector y.
B.3.2 Solution for S_{ik}(u), u ∈ U_{m}
Substituting c_{1,ik}, c_{2,ik}, c_{3,ik} from Equation A7 in Equation A4 we obtain for S_{ik}(u) :
B.3.3 Solution for P_{ik}(u), u ∈ U_{m}
Substituting c_{1,ik}, c_{2,ik}, c_{3,ik} from Equation A7 in Equation A4 we obtain for P_{ik}(u) :
B.3.4 Solution for D_{ik}(u), u ∈ U_{m}
Using Equation A1, Equation A8, and Equation A9 we obtain:
C Online supplement
C.1 Kaplan-Meier curves
C.2 JAGS code random effects NMA model
Example JAGS code corresponding to model SP Weibull RE2; SD exponential; PD Weibull FE1(scale) has been presented here.
model{ # Likelihood for conditional survival probabilities for (i in 1:Nd){ r_cond_pfs1[i]∼dbinom(cond_pfs[i,1], n_cond_pfs1[i]) r_cond_pfs2[i]∼dbinom(cond_pfs[i,2], n_cond_pfs2[i]) r_cond_os1[i]∼dbinom(cond_os[i,1], n_cond_os1[i]) r_cond_os2[i]∼dbinom(cond_os[i,2], n_cond_os2[i]) } # transformation of conditional survival probabilities to survival probabilities p4[1]<-1 p5[1]<-0 for (i in 2:Nd){ p4[i]<-equals(a[i]-a[i-1],0)*p[(i-1),4] + (1-equals(a[i]-a[i-1],0))*1 p5[i]<-equals(a[i]-a[i-1],0)*p[(i-1),5] + (1-equals(a[i]-a[i-1],0))*0 } for (i in 1:Nd){ cond_pfs[i,1]<-p[i,1]/p4[i] cond_pfs[i,2]<-p[i,4]/p4[i] cond_os[i,1]<-(p[i,1]+p[i,2])/(p4[i]+p5[i]) cond_os[i,2]<-(p[i,4]+p[i,5])/(p4[i]+p5[i]) # transformation of survival probabilities in time-varying hazards # for transitions between health states p[i,1]<-p4[i]*exp(-(h.sd[i]+h.sp[i])*dt[i,1]) p[i,2]<- p5[i]*exp(-h.pd[i]*dt[i,1])+p4[i]*h.sp[i]*(exp(-(h.sd[i]+h.sp[i])*dt[i,1]) -exp(-h.pd[i]*dt[i,1]))/(h.pd[i]-h.sp[i]-h.sd[i]) p[i,3]<-1-(p[i,1]+p[i,2]) p[i,4]<- p4[i]*exp(-(h.sd[i]+h.sp[i])*dt[i,2]) p[i,5]<- p5[i]*exp(-h.pd[i]*dt[i,2])+p4[i]*h.sp[i]*(exp(-(h.sd[i]+h.sp[i])*dt[i,2]) -exp(-h.pd[i]*dt[i,2]))/(h.pd[i]-h.sp[i]-h.sd[i]) p[i,6]<-1-(p[i,4]+p[i,5]) # decribe hazards as a function of time log(h.sp[i])<- Beta[s[i],a[i],1]+Beta[s[i],a[i],2]*timetrans1[i] log(h.sd[i])<- Beta[s[i],a[i],3] log(h.pd[i])<- Beta[s[i],a[i],4]+Beta[s[i],a[i],5]*timetrans1[i] } # random effects model for (l in 1:Ns){ for (ll in 1:na[l]){ Beta[l,ll,1]<-mu[l,1]+delta[l,ll] Beta[l,ll,2]<-mu[l,2]+d[t[l,ll],2]-d[t[l,1],2] Beta[l,ll,3]<-mu[l,3] Beta[l,ll,4]<-mu[l,4]+d[t[l,ll],3]-d[t[l,1],3] Beta[l,ll,5]<-mu[l,5] } w[l,1]<-0 delta[l,1]<-0 for (ll in 2:na[l]){ delta[l,ll]∼dnorm(md[l,ll],taud[l,ll]) md[l,ll]<-d[t[l,ll],1]-d[t[l,1],1] +sw[l,ll] w[l,ll] <-(delta[l,ll] - d[t[l,ll],1] + d[t[l,1],1]) sw[l,ll] <-sum(w[l,1:(ll-1)])/(ll-1) taud[l,ll] <-tau *2*(ll-1)/ll } } # priors for (j in 1:Ns){ mu[j,1:5] ∼ dmnorm(prior_mean_mu[1:5],prior_varcov_mu[,]) } d[1,1]<-0 d[1,2]<-0 d[1,3]<-0 for (k in 2:Nt){ d[k,1:3] ∼ dmnorm(prior_mean_d[1:3],prior_varcov_d[,]) } sd∼dunif(0,2) tau<-1/(sd*sd) # output for (k in 2:Nt){ for (u in 1:48){ log(HR.SP[1,k,u])<-(d[k,1]-d[1,1])+(d[k,2]-d[1,2])*log(u) log(HR.SD[1,k,u])<-0 log(HR.PD[1,k,u])<-(d[k,3]-d[1,3]) } } mu_mean[1]<--3.515819945 mu_mean[2]<-0.483167097 mu_mean[3]<--5.629762593 mu_mean[4]<--2.920038722 mu_mean[5]<--0.036010675 for (k in 1:Nt){ beta1[k]<-mu_mean[1]+d[k,1] beta2[k]<-mu_mean[2]+d[k,2] beta3[k]<-mu_mean[3] beta4[k]<-mu_mean[4]+d[k,3] beta5[k]<-mu_mean[5] } for (k in 1:Nt){ for (u in 1:48){ log(HAZARD.SP[k,u])<-(beta1[k])+(beta2[k])*log(u) log(HAZARD.SD[k,u])<-(beta3[k]) log(HAZARD.PD[k,u])<-(beta4[k])+(beta5[k])*log(u) }} for (k in 1:Nt){ P.S[k,1]<-1*exp(-(HAZARD.SD[k,1]+HAZARD.SP[k,1])) P.P[k,1]<-0*exp(-HAZARD.PD[k,1])+1*HAZARD.SP[k,1]*(exp(-(HAZARD.SD[k,1]+HAZARD.SP[k,1])) -exp(-HAZARD.PD[k,1]))/(HAZARD.PD[k,1]-HAZARD.SP[k,1]-HAZARD.SD[k,1]) PFS[k,1]<-P.S[k,1] OS[k,1]<-P.S[k,1]+P.P[k,1] for (u in 2:48){ P.S[k,u]<-P.S[k,(u-1)]*exp(-(HAZARD.SD[k,u]+HAZARD.SP[k,u])) P.P[k,u]<-P.P[k,(u-1)]*exp(-HAZARD.PD[k,u])+P.S[k,(u-1)]*HAZARD.SP[k,u] *(exp(-(HAZARD.SD[k,u]+HAZARD.SP[k,u])) -exp(-HAZARD.PD[k,u]))/(HAZARD.PD[k,u]-HAZARD.SP[k,u]-HAZARD.SD[k,u]) PFS[k,u]<-P.S[k,u] OS[k,u]<-P.S[k,u]+P.P[k,u] } } }