Bias and precision of methods for estimating the difference in restricted mean survival time from an individual patient data meta-analysis

Background The difference in restricted mean survival time (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗), the area between two survival curves up to time horizon \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {t}^{\ast } $$\end{document}t∗, is often used in cost-effectiveness analyses to estimate the treatment effect in randomized controlled trials. A challenge in individual patient data (IPD) meta-analyses is to account for the trial effect. We aimed at comparing different methods to estimate the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗ from an IPD meta-analysis. Methods We compared four methods: the area between Kaplan-Meier curves (experimental vs. control arm) ignoring the trial effect (Naïve Kaplan-Meier); the area between Peto curves computed at quintiles of event times (Peto-quintile); the weighted average of the areas between either trial-specific Kaplan-Meier curves (Pooled Kaplan-Meier) or trial-specific exponential curves (Pooled Exponential). In a simulation study, we varied the between-trial heterogeneity for the baseline hazard and for the treatment effect (possibly correlated), the overall treatment effect, the time horizon \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {t}^{\ast } $$\end{document}t∗, the number of trials and of patients, the use of fixed or DerSimonian-Laird random effects model, and the proportionality of hazards. We compared the methods in terms of bias, empirical and average standard errors. We used IPD from the Meta-Analysis of Chemotherapy in Nasopharynx Carcinoma (MAC-NPC) and its updated version MAC-NPC2 for illustration that included respectively 1,975 and 5,028 patients in 11 and 23 comparisons. Results The Naïve Kaplan-Meier method was unbiased, whereas the Pooled Exponential and, to a much lesser extent, the Pooled Kaplan-Meier methods showed a bias with non-proportional hazards. The Peto-quintile method underestimated the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗, except with non-proportional hazards at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {t}^{\ast } $$\end{document}t∗= 5 years. In the presence of treatment effect heterogeneity, all methods except the Pooled Kaplan-Meier and the Pooled Exponential with DerSimonian-Laird random effects underestimated the standard error of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗. Overall, the Pooled Kaplan-Meier method with DerSimonian-Laird random effects formed the best compromise in terms of bias and variance. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast },=,10,\kern0.5em ,\mathrm{years}\right) $$\end{document}rmstDt∗=10years estimated with the Pooled Kaplan-Meier method was 0.49 years (95 % CI: [−0.06;1.03], p = 0.08) when comparing radiotherapy plus chemotherapy vs. radiotherapy alone in the MAC-NPC and 0.59 years (95 % CI: [0.34;0.84], p < 0.0001) in the MAC-NPC2. Conclusions We recommend the Pooled Kaplan-Meier method with DerSimonian-Laird random effects to estimate the difference in restricted mean survival time from an individual-patient data meta-analysis. Electronic supplementary material The online version of this article (doi:10.1186/s12874-016-0137-z) contains supplementary material, which is available to authorized users.


Background
In cost-effectiveness analysis, a commonly used survival measure is the restricted mean survival time (RMST). It estimates the life expectancy for one treatment arm up to a certain time horizon t Ã [1][2][3][4]. The difference in restricted mean survival time (rmstD t Ã ð Þ) can thus quantify the treatment effect expressed in terms of life years gained. The rmstD t Ã ð Þ is an appealing outcome measure as it is valid even in case of non-proportional hazards [5]. Moreover, as an absolute outcome, the interpretation of the rmstD t Ã ð Þ is considered more intuitive from the clinician point of view than the hazard ratio (HR) which is a relative measure of the treatment effect [3,5,6]. Recent studies have compared methods to estimate the RMST including extrapolation beyond the trial followup [7][8][9][10]. However, these studies focused on the use of a single randomized clinical trial and not specifically on multicenter clinical trials nor meta-analyses.
In meta-analyses or in multicenter clinical trials, there is a need to take into account the trial or center effect to avoid the Simpson's paradox that may lead to biased estimates [11][12][13]. Different authors have discussed the use of Cox models with either stratification or fixed effect, or random effects to account for the center effect in a multicenter clinical trial [14][15][16] or the trial effect in a meta-analysis [17][18][19][20] in presence of baseline hazard heterogeneity and/or treatment effect heterogeneity between centers or trials. Several papers have also compared one-stage or two-stage methods to estimate the hazard ratio in individual patient data (IPD) metaanalyses [20][21][22]. All these studies focused on the estimation of the treatment effect through the use of the hazard ratio, but so far only one has focused on the use of the rmstD t Ã ð Þ in IPD meta-analyses [6]. In this latter study, Wei and colleagues investigated three two-stage methods to estimate the rmstD t Ã ð Þ from an IPD metaanalysis: two non-parametric methodsone based on pseudo-values [23] and one based on the Kaplan-Meier estimateand a flexible parametric survival model [24]. In their study, the rmstD t Ã ð Þ was estimated as an aggregation of the rmstD t Ã ð Þ estimated in each trial using a fixed effect meta-analysis model. The authors showed via simulations and two case studies that the three methods produced similar results in terms of bias and coverage probability of the confidence intervals.
In the present paper, we aimed at extending the first study from Wei et al. [6] by comparing more methods for estimating the rmstD t Ã ð Þ from an IPD meta-analysis in a broader range of scenarios. We also designed a more realistic simulation study with random effects to induce between-trial heterogeneity, both in terms of baseline hazard and of treatment effect. We considered only one of the non-parametric methods studied by Wei and colleaguesthe one pooling Kaplan-Meier estimates -as they found similar results for the three methods. We further considered a parametric methodpooling exponential estimatesand two other non-parametric methods: one naïve method that does not account for trial effect and an actuarial survival method developed by Peto, classically used in IPD meta-analyses for computing survival curves [25][26][27]. In simulations, we varied not only the treatment effect size and the time horizon t Ã , as previously done by Wei and colleagues, but also the baseline hazard heterogeneity, the treatment effect heterogeneity, the correlation between these two random effects, the number of trials and the number of patients per trial, the use of fixed effect and DerSimonian-Laird random effects model [28], and departure from the assumption of proportional hazards.
In the 'Methods' section we describe the rmstD t Ã ð Þ and the different survival analysis methods for estimating the rmstD t Ã ð Þ that we investigate in this paper. Section 'Simulation study' describes the design of the simulation study, how to estimate the true rmstD t Ã ð Þ, the simulation scenarios and the evaluation criteria, and presents the simulation results. Section ' Application' provides two examples using IPD meta-analyses in nasopharynx carcinoma. We end with a discussion and our conclusion regarding the choice of a particular method to estimate the rmstD t Ã ð Þ from an IPD meta-analysis. Of note, the investigated methods can also be used for the estimation of the rmstD t Ã ð Þ in multicenter clinical trials.

Difference in restricted mean survival time
Let T be the survival time random variable with distribution F(t). The mean survival time restricted at a specified time t Ã is defined as where S(t) = 1 -F(t) is the survival function. The RMST t Ã corresponds graphically to the area under the survival curve S(t) up to t Ã . The difference in restricted mean survival time (rmstD t Ã ð Þ) between the experimental arm and the control arm (noted 1 and 0) is defined as The variance d V ar rmstD t Ã ð Þ ð Þ can be estimated as [29]: As opposed to the relative hazard ratio, the rmstD t Ã ð Þ is an absolute outcome which depends both on the baseline hazard and on the relative treatment effect.
Wei et al. [6] also proposed the use of the relative difference in restricted mean survival time defined as The rmstD t Ã ð Þ varies between 0 and 1 and can be interpreted as a percentage. Its variance can be esti- Methods for estimating the difference in restricted mean survival time We investigated four methods for estimating the difference in restricted mean survival time rmstD t Ã ð Þ from an IPD meta-analysis: 1) the Naïve Kaplan-Meier, which pools the data, ignoring the trial effect, 2-3) the Pooled Kaplan-Meier and Pooled Exponential methods, which use a two-stage approach to combine rmstD j t Ã ð Þ estimated in each trial j, and 4) the Peto-quintile method, which uses survival functions derived from a pooled hazard ratio estimated with a two-stage approach in order to take into account the trial effect.

Naïve Kaplan-Meier
The most basic method to estimate the rmstD t Ã ð Þ is to consider the IPD meta-analysis as a single large trial. Under this approach, the rmstD t Ã ð Þ is estimated as the area between the Kaplan-Meier curves of the experimental and the control arm, obtained by pooling the data from all the trials, thus ignoring the trial effect [30]: with Ŝ arm (t arm,0 ) = 0, t arm,0 = 0, D arm the number of distinct event times (t arm,1 < t arm,2 < … < t arm,D ) and Ŝ arm (t) the Kaplan-Meier estimator for the experimental arm and the control arm (noted 1 and 0). The variance of the rmstD t Ã ð Þ was estimated analytically by the delta method for the Naïve Kaplan-Meier (details are provided in Additional file 1).

Pooled Kaplan-Meier and Pooled Exponential
In order to take into account the trial effect, a different strategy consists in estimating the rmstD j t Ã ð Þ within each trial j and then to combine the trial-specific results into a pooled estimate. In the Pooled Kaplan-Meier and Pooled Exponential methods, which are both two-stage approaches, we first estimated the rmstD j t Ã ð Þ in each trial j as the area between either trial-specific Kaplan-Meier curves or trial-specific exponential curves with for each trial j: Ŝ j,arm (t j,arm,0 ) = 0, t j,arm,0 = 0, D j,arm the number of distinct event times (t j,arm,1 < t j,arm,2 < … < t j,arm,D ), Ŝ j,arm (t) the Kaplan-Meier estimator, andλ j;arm the maximum likelihood estimate of the scale parameter for the exponential distribution for the experimental arm and the control arm (noted 1 and 0). Then, we combined the rmstD j t Ã ð Þ by using a fixed effect or a DerSimonian-Laird random effects [28] meta-analysis model (see Additional file 1). The variance of each rmstD j t Ã ð Þ was estimated analytically by the delta method for the Pooled Kaplan-Meier and Pooled Exponential methods (details are provided in Additional file 1).

Peto-quintile
The actuarial method was developed by Peto to compute the survival curves taking into account the trial effect in IPD meta-analyses [25][26][27]31]. In this case, the survival probabilities in the two arms are estimated at the end of predetermined time intervals i based on the estimated survival probability p i in the overall population and a pooled hazard ratio HR i,pooled . The survival probability for the overall population is estimated as where D i is the number of deaths during interval i and PI i the total number of person-intervals in the i-th interval. One person-interval is equivalent to one person-year when the time interval chosen is 1 year. The pooled hazard ratio in the interval, HR b i;pooled ; is estimated using a fixed effect meta-analysis model to aggregate the different HR b i;j estimated in each trial j. The survival probabilities at the end of each interval i in the control arm (p 0,i ) and in the experimental arm (p 1,i ) are estimated as follows: The survival probability at time t in each arm is the product of the survival probabilities across n i intervals up to t In the present work, we extended this method by also pooling the d HRi; j using a DerSimonian-Laird random effects meta-analysis model (see Additional file 1). Furthermore, we used time intervals i based on the quintiles of the overall number of deaths occurring before t Ã in the meta-analysis and therefore we called this method the Peto-quintile method. The rmstD t Ã ð Þ was then estimated as the area between the two survival curves defined by Ŝ Peto,0 (t) and Ŝ Peto,1 (t).
where t 0 = 0 and t 1 ;::; Þdenotes the time intervals based on the quintiles of events.
A 50-replicate non-parametric bootstrap was used to estimate the variance of the rmstD t Ã ð Þ for the Petoquintile method.

Follow-up differences across trials and extrapolation
We used the extrapolation method proposed by Brown et al. [32] for Naïve Kaplan-Meier and Pooled Kaplan-Meier, to extrapolate the survival function beyond the last observed event time (t max ) until t Ã , if needed (e.g. in case of potential follow-up difference across trials) [10,32]. The Brown extrapolation method consists in completing the tail of the Kaplan-Meier survival curve by an exponential curve. The estimated survival function for t > t max is: where Ŝ (t) is the Kaplan-Meier estimator of the survival function.
No extrapolation was needed for Pooled Exponential as a parametric model is used. Concerning the Petoquintile method, at least one event is needed overall in the meta-analysis in each arm and at each time interval i for survival probabilities p 0,i and p 1,i to be computed. This was always the case, even in case of a potential difference in follow-up across trials, as each time interval contained by definition one fifth of the deaths occurring before t Ã . It is worth noting that if t Ã is greater than t max the last observed event in the whole meta-analysis, the method estimates the event rate p i and the pooled hazard ratio HR i,pooled until t max and implicitly extrapolates the survival in the remaining interval t max ; t Ã ½ .

Simulation study
Design of the simulation study We simulated survival times of N patients from J trials, each of size n j with ∑ j = 1 J n j = N. We defined the hazard function for a trial j ϵ {1 , …, J} as where λ 0 (t) is the baseline hazard function, A j is a trialspecific random quantity affecting the baseline hazard with Var(A j ) = σ 2 , and B j is a trial-specific random quantity affecting the treatment effect with Var(B j ) = τ 2 , β is the overall treatment effect, and χ is the binary treatment variable, coded +1/2 for experimental arm, −1/2 for control in order to obtain equal heterogeneity in both arms [19,21,22]. Of note, the covariance between the two random effects is defined by cov(A j ,B j ) = ρστ with ρ the correlation between A j and B j . We used an exponential distribution for the baseline hazard function λ 0 (t) = (log(2)/5)t, corresponding to a median survival time of 5 years. Independent and noninformative right-censoring was induced by setting the recruitment time at 3 years for all trials and varying the maximum follow-up time uniformly between 2 and 9 years across trials to replicate the typical difference in observed follow-up between trials included an IPD meta-analysis.
We induced between-trial heterogeneity by generating random values a j and b j from binomial distributions for the baseline hazard and the treatment effect. The use of a discrete distribution allowed us to derive straightforwardly the true difference in restricted mean survival time ( rmstD t Ã ð Þ). The binomial random variables were centered and properly rescaled in order to obtain the desired variances σ 2 and τ 2 : The rationale for the arbitrary choice of n = 50 was that the distribution approximated well a continuous distribution, while allowing easy computation of the true rmstD t Ã ð Þ.

True difference in restricted mean survival time
Based on our simulation model, the difference in restricted mean survival time is defined as Thanks to the use of a discrete joint probability distribution F A,B (a,b), the integral in equation (18) boils down to a sum over the K points belonging to its support K: with the conditional restricted mean survival time rmstD k t Ã ð Þ defined for a couple (a k, b k ) as: Simulation scenarios In different scenarios we varied the strength of betweentrial heterogeneity for the baseline hazard (low with σ 2 = 0.01 and high with σ 2 = 0.10) and for the treatment effect (τ 2 = 0.01, 0.10). We performed the main analysis with uncorrelated random effects (ρ = 0), however, in a sensitivity analysis, we studied the impact of a negative correlation between A j and B j (ρ = −0.8). We considered different values for the number of trials and patients per trial: (J = 5, n j = 200) and (J = 20, n j = 100) and for the size of the overall treatment effect (β = 0, ±0.2, ±0.7). We also studied the impact of the time horizon of restriction at t Ã = 5 years and t Ã = 10 years. These two values were chosen to illustrate a scenario in which all trials still have patients at risk at t Ã (5 years) and a scenario in which some trials' follow-up are shorter than the time of restriction ( t Ã = 10 years). The average administrative censoring rate ranged across scenarios from 49 to 52 % at t Ã = 5 years and from 38 to 40 % at t Ã = 10 years. In the case of no overall treatment effect at all (β = 0, τ 2 = 0) and no baseline heterogeneity (σ 2 = 0), the restricted mean survival time was equal to 3.6 years at t Ã = 5 years and 5.4 years at t Ã = 10 years in both arms. The influence of non-proportional hazards was examined using a piecewise exponential distribution with a deleterious treatment effect (β' = −β, with β ≤ 0) in the first 2 years and a beneficial treatment effect (β) afterwards.

Evaluation criteria
We simulated 1,000 meta-analyses for each scenario and compared the four methods using: the average bias, defined as the average of the estimated rmstD t Ã ð Þ minus the true value; the empirical standard error (ESE), defined as the standard deviation of the rmstD t Ã ð Þ over the replicates; and the average standard error (ASE), defined as the average of the estimated standard errors [33].
With the exception of the Naïve Kaplan-Meier method, the methods were not available in standard statistical software. We implemented the methods and performed the simulation study using R version 3.1.3 (R Foundation, Vienna, Austria). The R code is available from the authors upon request.

Results
For all scenarios, there was almost no bias in the case of no treatment effect (β = 0). When there was a beneficial treatment effect (β = −0.2 or −0.7), Peto-quintile underestimated the rmstD t Ã ð Þ, notably on the long term (t Ã = 10 years). The Pooled Exponential and, to a much lesser extent, the Pooled Kaplan-Meier methods showed a bias in the case of non-proportional hazards. In all these cases, whenever a method showed a bias, the bias increased with |β| (Table 1 and Fig. 1).
In scenarios with higher treatment effect heterogeneity (τ 2 = 0.10), all the methods had higher empirical standard error (ESE), as shown in the Figs. 1 and 2, and Tables 1 and 2. The standard error was estimated correctly (ASE = ESE) only with Pooled Kaplan-Meier and Pooled Exponential. It was generally underestimated (ASE < ESE) with the two other methods: the ASE was two-fold smaller than the ESE for the Naïve Kaplan-Meier and the Peto-quintile methods with τ 2 = 0.10. When varying the baseline hazard heterogeneity between trials, no relevant impact was noted neither on the bias nor on the standard error.
With both proportional and non-proportional hazards, for β = −0.7, the Peto-quintile method showed a bias which was negligible at t Ã = 5 years but much higher at t Ã = 10 years (up to 0.21 years; Fig. 1-b and Fig. 2-b). In the case of non-proportional hazards, which were incorporated using a piecewise exponential distribution with a deleterious treatment effect in the first 2 years and a beneficial treatment effect afterwards, Pooled Exponential was heavily biased at 5 years, with a bias of almost 0.40 years as compared to a true rmstD t Ã ¼ 5 ð Þ= −0.30 years (Fig. 2-a and Table 2). This bias suggests that the Pooled Exponential method failed to reflect the piecewise exponential distribution with β' = −β (β ≤ 0) for t ϵ [0;2] years and β for t > 2 years. However, this bias disappeared at 10 years, arguably because the true hazards were proportional between 2 and 10 years in our simulation set-up ( Fig. 2-b) and the different effect in the first 2 years was thus attenuated. A small bias also arose for Pooled Kaplan-Meier at 10 years when the hazards were not proportional: a bias of around 0.05 years as compared to a true rmstD t Ã ¼ 10 ð Þ= 0.30 years ( Fig. 2-b). In terms of standard error, lower values were found for both ESE and ASE at t Ã = 5 years ( Fig. 1-a and Fig. 2-a) as compared at t Ã = 10 years ( Fig. 1-b and Fig. 2-b), no matter the hazards were proportional or not.
The number of trials and the size of trials had no major impact in terms of bias. In meta-analyses of J = 20 trials and n j = 100 patients per trial, all the methods had lower empirical and average standard errors than in meta-analyses of J = 5 trials and n j = 200 patients per trial (see Additional file 1: Tables S1, 2, 3 and 4).
We also considered a deleterious treatment effect (β = +0.2, +0.7) but, as expected, results were not affected: the biased methods had biases that were reversed, and ASE and ESE remained unchanged (Additional file 1: Table S5).
The introduction of a negative correlation between A j and B j also had no major impact in terms of bias and standard error estimation, with the exception of the scenario with high baseline hazard and treatment variances (σ 2 = τ 2 = 0.10) for which ASE and ESE of all the methods were higher than with no correlation, notably for β = −0.7 (Additional file 1: Table S6).
When using a fixed effect meta-analysis model (Additional file 1: Table S7), for scenarios with high treatment effect heterogeneity (τ 2 = 0.10), the Pooled Kaplan-Meier, Pooled Exponential and Peto-quintile methods exhibited a larger bias as compared with a DerSimonian-Laird random effects models used in Table 1. Furthermore, using a fixed effect model underestimated the standard error in general (ASE < <ESE).

Application
We illustrate the four methods for estimating the rmstD t Ã ð Þ using IPD from the Meta-Analysis of Chemotherapy in Nasopharynx Carcinoma (MAC-NPC) Collaborative Group [34] and its updated version MAC-NPC2 [35] as these two IPD meta-analyses differed in terms of evidence of treatment effect heterogeneity. These IPD meta-analyses studied the addition of chemotherapy (CT) to radiotherapy (RT) in patients with nasopharynx carcinoma. For the estimation of the rmstD t Ã ð Þ, we selected t Ã = 5 years and t Ã = 10 years, as these were the two time points of clinical interest in the publications of MAC-NPC and MAC-NPC2. The overall proportional hazards assumption was verified at the 5 % significance level (p = 0.09) according to the methodology described by Wei et al. [6], in which trial-specific p-values from Grambsch-Therneau test [36] are pooled. The rmstD t Ã ð Þ ranged from 0.17 to 0.23 years at t Ã = 5 years and from 0.46 to 0.55 years at t Ã = 10 years across the estimation methods (Table 3). For Pooled Kaplan-Meier and Pooled Exponential using a random effects model, the rmstD t Ã ð Þ was not significantly different from 0. As there was high treatment effect heterogeneity in the MAC-NPC, a DerSimonian-Laird random effects model was deemed more appropriate to aggregate the trial-specific rmstD j t Ã ð Þ. As previously seen in the simulation study, a fixed effect model would underestimate the variance of the overall estimate. Also, similarly to our simulation study with proportional hazards, larger values for rmstD t Ã ð Þ and SE (rmstD t Ã ð Þ) were found at t Ã = 5 years as compared to at t Ã = 10 years. Figure 3 displays the forest plot for trial-specific rmstD j t Ã ð Þ and overall rmstD t Ã ð Þ estimated using Pooled Kaplan-Meier with DerSimonian-Laird random effects at t Ã = 10 for the MAC-NPC metaanalysis. Figure 4 displays the overall rmstD(t * ) estimated by Pooled Kaplan-Meier with DerSimonian-Laird random effects when varying t Ã ; it shows that the rmstD t Ã ð Þ is not significantly different from 0 for t Ã ϵ [0;10] years. The same graphic for the overall rmstD t Ã ð Þ is displayed as Additional file 1: Figure S1.  The pooled p-value test (p = 0.16) suggested that the overall proportional hazards assumption was appropriate. The rmstD t Ã ð Þ ranged from 0.17 to 0.21 years at t Ã = 5 years and from 0.53 to 0.59 years at t Ã = 10 across the estimation methods ( Table 3). The rmstD t Ã ð Þ was significantly different from 0 and in favor of the RT + CT arm with all of the methods. As compared to the results in the MAC-NPC, the standard error of the rmstD t Ã ð Þ was lower in the MAC-NPC2 with a rmstD t Ã ð Þ of similar magnitude for all the methods. This was consistent with the simulation results, as there were more trials and overall more patients included in the MAC-NPC2. The forest plot for the MAC-NPC2 displaying trial-specific rmstD j t Ã ð Þ and the overall rmstD t Ã ð Þ at t Ã = 10 years estimated using the Pooled Kaplan-Meier method with DerSimonian-Laird random effects is provided in Additional file 1: Figure S2.

Discussion
The difference in restricted mean survival time (rmstD t Ã ð Þ) is an appealing alternative to the hazard ratio (HR) as measure of treatment effect, because it does not require the proportional hazards assumption and is considered to have a  MAC-NPC meta-analysis of chemotherapy in nasopharynx carcinoma, rmstD difference in restricted mean survival time, SE standard error, t * time horizon more intuitive interpretation [3,5,6]. Furthermore, the rmstD t Ã ð Þ is directly related to cost-effectiveness analysis as it is the denominator of the incremental cost-effectiveness ratio, so one can use the rmstD t Ã ð Þ estimation from a previous publication to perform a cost-effectiveness analysis. We previously showed that in a cost-effectiveness analysis even small variations in the estimate of the rmstD t Ã ð Þ from an individual patient data (IPD) meta-analysis can yield significantly different reimbursement conclusions [37]. However, to our knowledge only one evaluation of the methods to estimate the rmstD t Ã ð Þ from IPD meta-analysis is available to date [6].
In this study, we compared different methods to estimate the rmstD t Ã ð Þ from IPD meta-analysis in different scenarios varying several key meta-analysis parameters. We showed that Pooled Kaplan-Meier was rarely biased. Similarly, Naïve Kaplan-Meier was unbiased in all scenarios, whereas Pooled Exponential showed a bias with non-proportional hazards at t Ã = 10 years and an even larger bias at t * = 5 years. Peto-quintile underestimated the rmstD t Ã ð Þ , except with non-proportional hazards at t Ã = 5 years. In case of treatment effect heterogeneity, the use of a fixed effect model was not appropriate and all methods except Pooled Kaplan-Meier and Pooled Exponential with DerSimonian-Laird The rmstD t Ã ð Þ is an absolute outcome measure which depends both on the baseline hazard and on the relative treatment effect. Consequently, the heterogeneity test when pooling the rmstD j t Ã ð Þ reflects both baseline hazard and relative treatment effect heterogeneities. Deeks already showed that in 551 systematic reviews with binary outcomes the heterogeneity was higher for an absolute outcome than for a relative outcome [38]. In our IPD metaanalyses in nasopharynx carcinoma, the heterogeneity was slightly higher when pooling the rmstD j t Ã ð Þ with Pooled Kaplan-Meier than when pooling the hazard ratios: for the MAC-NPC there was a small increase in the heterogeneity with Cochran Q test p-value = 0.03, I 2 = 50 % for HR as compared to p = 0.02, I 2 = 54 % for rmstD t Ã ¼ 10 years ð Þ (Fig. 3). For the MAC-NPC2, this increase was more pronounced with p = 0.09, I 2 = 30 % for HR versus p = 0.01, I 2 = 45 % for rmstD t Ã ¼ 10 years ð Þ(Additional file 1: Figure  S1). Wei and colleagues showed a similar trend in their second example (p = 0.47, I 2 = 0 % for HR and p = 0.20, I 2 = 24 % for rmstD t Ã ¼ 5 years ð Þ ) [6]. In our simulation study, we have induced between-trial heterogeneity for the baseline hazard and for the treatment effect using two random effects. As a matter of fact, both random effects can be tested, by testing Var(A j ) = σ 2 = 0 and Var(B j ) = τ 2 = 0. Testing for the presence of treatment effect heterogeneity (τ 2 = 0) corresponds to the Cochran Q-test which we have used in the MAC-NPC applications [34,35]. Commenges and Andersen [39] and Biard and colleagues [40] proposed respectively the use of score tests or permutation tests for testing the baseline heterogeneity (σ 2 = 0) in proportional hazard models. Rondeau et al. [19] tested both the baseline hazard (σ 2 = 0) and the treatment effect heterogeneity between trials (τ 2 = 0) using a mixture of χ 2 distributions in one-stage Cox models.
Recent techniques like the one proposed by Guyot et al. [41] allow one to reconstruct IPD based on published Kaplan-Meier curves, which could be useful to recalculate the rmstD t Ã ð Þ even for aggregate data. However, we suggest that clinical publications for single (multicenter) clinical trial or IPD meta-analysis should report the rmstD t Ã ð Þ at different time horizons t Ã of clinical interest in addition to the hazard ratio. This way the rmstD t Ã ð Þ would be available for future economic evaluations. This is of particular relevance as a previous study stated that the survival outcome in a cost-effectiveness analysis based on a clinical trial or a meta-analysis should be estimated with the same statistical model used for efficacy [42].
Among the two-stage methods studied by Wei et al., the non-parametric pseudo-values method was disregarded, as Wei et al. showed that it led to similar results as the nonparametric Pooled Kaplan-Meier method [6]. Also, among parametric models we chose the exponential model instead of the Royston and Parmar flexible parametric model for ease of computation. In addition, we chose to study other non-parametric methods from the medical literature that have been actually applied in practice. Parametric methods developed for network meta-analysis were not included [43,44]. Furthermore, methods using the percentile ratio [45,46] were beyond the scope of this study, which focused on the rmstD t Ã ð Þ. In addition, in this simulation study, we only considered balanced trials and we did not vary the administrative censoring rate.
The rmstD t Ã ð Þ is inherently dependent on the choice of t Ã . Also, we showed that its standard error gets larger as t Ã increases (Fig. 1). Karrison recommended to choose a maximum time horizon t Ã such that SE(S( t Ã )) is less than a chosen ceiling value [29,47]. In the particular case of an IPD meta-analysis, trials can have different lengths of follow-up, and there is thus a compromise to achieve between small values of t Ã that censor a lot of data with a high loss of information, and high values of t Ã that need a massive use of extrapolation. Wei and colleagues stated that the choice of t Ã should also be of clinical interest, and they suggested plotting the rmstD t Ã ð Þ against t Ã to see how the treatment effect varies over time. In MAC-NPC for instance such a plot shows that the rmstD t Ã ð Þ was not significantly different from 0 with t Ã ϵ [0;10] years based on pointwise confidence intervals (Fig. 4). In two recent papers, Tian et al. [48] and Zhao et al. [4] have proposed a simultaneous confidence interval of the rmstD t Ã ð Þ in the context of one randomized controlled trial. However, an extension to the context of IPD meta-analyses or multicenter clinical trials has not yet been proposed and may be the subject of further research.
Depending on the choice of the time horizon t Ã , some trials included in the IPD meta-analysis may have a follow-up not long enough to reach t Ã . In our study, for such trials, we used the extrapolation method proposed by Brown et al. [32] until t Ã for the Naïve Kaplan-Meier and the Pooled Kaplan-Meier methods. Lamb and colleagues [10] have shown that this extrapolation method is less biased than the mean survival time restricted at the last observed event time. For lifetime extrapolation, which can be needed in cost-effectiveness analysis, one can estimate the difference in mean survival time using the Pooled Kaplan-Meier with a DerSimonian-Laird random effects model. In each trial, the difference in mean survival time would be estimated using Kaplan-Meier curves with extrapolated parametric tails [9,10]. Similarly, for the two other non-parametric methods Naïve Kaplan-Meier and Peto-quintile, one can extrapolate the survival curves beyond the last observed failure time by using an extrapolated parametric tail.

Conclusions
The difference in restricted mean survival time (rmstD t Ã ð Þ) is an appealing alternative to the hazard ratio to measure the treatment effect in a meta-analysis of time-to-event outcomes, as it is free of the proportional hazards assumption and its interpretation is more intuitive. We compared methods to estimate the rmstD t Ã ð Þ from an individual patient data meta-analysis. In our simulation study, in which a large panel of meta-analysis parameters was varied, the two-stage Pooled Kaplan-Meier method with DerSimonian-Laird random effects formed the best compromise in terms of bias and variance. Thus, Pooled Kaplan-Meier with DerSimonian-Laird random effects should be the preferred method to estimate the difference in restricted mean survival time from an individual patient data meta-analysis or from a multicenter clinical trial.

Availability of data and materials
Data were used with permission obtained from the MAC-NPC Collaborative Group investigators, who agreed to share their data with us by signing an amendment to the original protocol. The French data protection authority (CNIL -Commission Nationale de l'Informatique et des Libertés) does not allow us to make these data publicly available.

Additional file
Additional file 1: Supplementary statistical details, Tables S1-S7 and

Competing interest
The authors declare that they have no competing interests.
Authors' contributions BL and FR performed the statistical analyses and drafted the manuscript. BL, FR, JB, JPP, and SM jointly contributed to writing the study protocol, interpreting and discussing results, and to writing the article. BL, FR, JB, JPP, and SM read and approved the final manuscript.