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

Abstract

Background

The difference in restricted mean survival time ($rmstD\left({t}^{\ast}\right)$), the area between two survival curves up to time horizon ${t}^{\ast }$, 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 $rmstD\left({t}^{\ast}\right)$ 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 ${t}^{\ast }$, 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 $rmstD\left({t}^{\ast}\right)$, except with non-proportional hazards at ${t}^{\ast }$= 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 $rmstD\left({t}^{\ast}\right)$. Overall, the Pooled Kaplan-Meier method with DerSimonian-Laird random effects formed the best compromise in terms of bias and variance. The $rmstD\left({t}^{\ast },=,10,\kern0.5em ,\mathrm{years}\right)$ 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.

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}^{\ast }$ [1–4]. The difference in restricted mean survival time ($rmstD\left({t}^{\ast}\right)$) can thus quantify the treatment effect expressed in terms of life years gained. The $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ 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 follow-up [7–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–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–16] or the trial effect in a meta-analysis [17–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) meta-analyses [20–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\left({t}^{\ast}\right)$ in IPD meta-analyses [6]. In this latter study, Wei and colleagues investigated three two-stage methods to estimate the $rmstD\left({t}^{\ast}\right)$ from an IPD meta-analysis: two non-parametric methods – one based on pseudo-values [23] and one based on the Kaplan-Meier estimate – and a flexible parametric survival model [24]. In their study, the $rmstD\left({t}^{\ast}\right)$ was estimated as an aggregation of the $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ 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 colleagues – the one pooling Kaplan-Meier estimates – as they found similar results for the three methods. We further considered a parametric method – pooling exponential estimates – and 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–27]. In simulations, we varied not only the treatment effect size and the time horizon ${t}^{\ast }$, 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\left({t}^{\ast}\right)$ and the different survival analysis methods for estimating the $rmstD\left({t}^{\ast}\right)$ that we investigate in this paper. Section ‘Simulation study’ describes the design of the simulation study, how to estimate the true $rmstD\left({t}^{\ast}\right)$, 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\left({t}^{\ast}\right)$ from an IPD meta-analysis. Of note, the investigated methods can also be used for the estimation of the $rmstD\left({t}^{\ast}\right)$ in multicenter clinical trials.

Methods

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}^{\ast }$ is defined as

$$ RMST\left({t}^{\ast}\right)={\displaystyle {\int}_0^{t^{\ast }}S(t)dt,} $$

(1)

where S(t) = 1 – F(t) is the survival function. The RMST ${t}^{\ast }$ corresponds graphically to the area under the survival curve S(t) up to ${t}^{\ast }$. The difference in restricted mean survival time ($rmstD\left({t}^{\ast}\right)$) between the experimental arm and the control arm (noted 1 and 0) is defined as

$$ rmstD\left({t}^{\ast}\right)=RMS{T}_1\left({t}^{\ast}\right)-RMS{T}_0\left({t}^{\ast}\right) $$

(2)

The variance $\widehat{Var}\left( rmstD\left({t}^{\ast}\right)\right)$ can be estimated as [29]:

$$ \widehat{Var}\left( rmstD\left({t}^{\ast}\right)\right)=\widehat{Var}\left(RMS{T}_1\left({t}^{\ast}\right)\right)+\widehat{Var}\left(RMS{T}_0\left({t}^{\ast}\right)\right) $$

(3)

As opposed to the relative hazard ratio, the $rmstD\left({t}^{\ast}\right)$ 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

$$ rmstRD\left({t}^{\ast}\right)= rmstD\left({t}^{\ast}\right)/{t}^{\ast } $$

(4)

The $rmstD\left({t}^{\ast}\right)$ varies between 0 and 1 and can be interpreted as a percentage. Its variance can be estimated as $\widehat{Var}\left[ rmstD\left({t}^{\ast}\right)\right]/{\left({t}^{\ast}\right)}^2$.

Methods for estimating the difference in restricted mean survival time

We investigated four methods for estimating the difference in restricted mean survival time $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ is to consider the IPD meta-analysis as a single large trial. Under this approach, the $rmstD\left({t}^{\ast}\right)$ 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]:

$$ \widehat{rmst}D\left({t}^{\ast}\right)={\displaystyle {\sum}_{i=1}^{D_1-1}{\widehat{S}}_1\left({t}_{1,i}\right)}\left({t}_{1,i+1}-{t}_{1,i}\right)-{{\displaystyle {\sum}_{i=1}^{D_0-1}\widehat{S}}}_0\left({t}_{0,i}\right)\left({t}_{0,i+1}-{t}_{0,i}\right) $$

(5)

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\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ in each trial j as the area between either trial-specific Kaplan-Meier curves

$$ \widehat{rmst}{D}_j\left({t}^{\ast}\right)={\displaystyle {\sum}_{i=1}^{D_{j,1}-1}{\widehat{S}}_{j,1}\left({t}_{j,1,i}\right)\left({t}_{j,1,i+1}-{t}_{j,1,i}\right)-{\displaystyle {\sum}_{i=1}^{D_{j,0}-1}{\widehat{S}}_{j,0}\left({t}_{j,0,i}\right)\left({t}_{j,0,i+1}-{t}_{j,0,i}\right)}} $$

(6)

or trial-specific exponential curves

$$ \widehat{rmst}{D}_j\left({t}^{\ast}\right)={\displaystyle {\int}_0^{t^{\ast }}{e}^{-{\widehat{\lambda}}_{j,1}t}dt-{\displaystyle {\int}_0^{t^{\ast }}{e}^{-{\widehat{\lambda}}_{j,0}t}dt=\frac{1-{e}^{-{\widehat{\lambda}}_{j,1}{t}^{\ast }}}{{\widehat{\lambda}}_{j,1}}}}-\frac{1-{e}^{-{\widehat{\lambda}}_{j,0}{t}^{\ast }}}{{\widehat{\lambda}}_{j,0}} $$

(7)

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 ${\widehat{\lambda}}_{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\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ 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–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

$$ \widehat{p_i}= \exp \left(-\frac{D_i}{P{I}_i}\right) $$

(8)

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, ${\widehat{HR}}_{i, pooled},$ is estimated using a fixed effect meta-analysis model to aggregate the different ${\widehat{HR}}_{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:

$$ \widehat{p_{0,i}}=\widehat{p_i}-\left[0.5\widehat{p_i}\left(\widehat{p_i}-1\right) log\left({\widehat{HR}}_{i,\kern0.1em pooled}\right)\right] $$

(9)

$$ \widehat{p_{1,i}}=\widehat{p_i}+\left[0.5\widehat{p_i}\left(\widehat{p_i}-1\right) log\left({\widehat{HR}}_{i,\kern0.1em pooled}\right)\right] $$

(10)

The survival probability at time t in each arm is the product of the survival probabilities across n _i intervals up to t

$$ {\widehat{S}}_{peto,0}(t)={\displaystyle {\prod}_{i=1}^{n_i}\hat{p_{0,i}}}\kern0.5em \mathrm{and}\kern0.5em {\widehat{S}}_{peto,1}(t)={\displaystyle {\prod}_{i=1}^{n_i}\hat{p_{1,i}}} $$

(11)

In the present work, we extended this method by also pooling the ${\widehat{HR}}_{i,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}^{\ast }$ in the meta-analysis and therefore we called this method the Peto-quintile method. The $rmstD\left({t}^{\ast}\right)$ was then estimated as the area between the two survival curves defined by Ŝ _Peto,0(t) and Ŝ _Peto,1(t).

$$ \widehat{rmstD}\left({t}^{\ast}\right)={\displaystyle {\sum}_{i=0}^4\frac{\left({t}_{i+1}-{t}_i\right)}{2}\left[\left({\widehat{S}}_{Peto,1}\left({t}_{i+1}\right)-{\widehat{S}}_{Peto,0}\left({t}_{i+1}\right)\right)+\left({\widehat{S}}_{Peto,1}\left({t}_i\right)-{\widehat{S}}_{Peto,0}\left({t}_i\right)\right)\right]} $$

(12)

where t ₀  = 0 and $\left({t}_1,..,\kern0.2em {t}_5={t}^{*}\right)$ 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\left({t}^{\ast}\right)$ for the Peto-quintile 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}^{\ast }$, 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:

$$ {\widehat{S}}_{Brown}(t)= \exp \left\{t\times log\left[\widehat{S}\left({t}_{max}\right)\right]/{t}_{max}\right\} $$

(13)

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 Peto-quintile 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}^{\ast }$. It is worth noting that if ${t}^{\ast }$ 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 $\left[{t}_{max},{t}^{*}\right]$.

Simulation study

Design of the simulation study

We simulated survival times of N patients from J trials, each of size n _j with ∑ ^J _j = 1 n _j  = N. We defined the hazard function for a trial j ϵ {1 , …, J} as

$$ \lambda \left(t;x\Big|{A}_j={a}_j,{B}_j={b}_j\right)={\lambda}_0(t) exp\left\{{a}_j+\left(\beta +{b}_j\right)x\right\} $$

(14)

where λ ₀ (t) is the baseline hazard function, A _j is a trial-specific random quantity affecting the baseline hazard with Var(A _j ) = σ ², and B _j is a trial-specific random quantity affecting the treatment effect with Var(B _j ) = τ ², β is the overall treatment effect, and x 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 λ₀(t) = (log(2)/5)t, corresponding to a median survival time of 5 years. Independent and non-informative 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\left({t}^{\ast}\right)$). The binomial random variables were centered and properly rescaled in order to obtain the desired variances σ ² and τ ²:

$$ A\sim \left[ Bin\left(n=50,p=0.5\right)-25\right]\cdot \sigma /\sqrt{12.5} $$

(15)

$$ B\sim \left[ Bin\left(n=50,p=0.5\right)-25\right]\cdot \tau /\sqrt{12.5} $$

(16)

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\left({t}^{\ast}\right)$.

True difference in restricted mean survival time

Based on our simulation model, the difference in restricted mean survival time is defined as

$$ rmstD\left({t}^{*}\right)={\displaystyle {\int}_0^{t^{*}}S\left(t;x=\frac{1}{2}\right)dt-}{\displaystyle {\int}_0^{t^{*}}S\left(t;x=-\frac{1}{2}\right)(t)dt} $$

(17)

$$ rmstD\left({t}^{*}\right)={\displaystyle {\int}_0^{t^{*}}\left[{\displaystyle {\int}_{\mathcal{K}}S}\left(t;x=\frac{1}{2}\Big|{a}_j,{b}_j\right)d{F}_{A,B}\left(a,b\right)\right]dt-{\displaystyle {\int}_0^{t^{*}}\left[{\displaystyle {\int}_{\mathcal{K}}S}\left(t;x=-\frac{1}{2}\Big|{a}_j,{b}_j\right)d{F}_{A,B}\left(a,b\right)\right]dt}} $$

(18)

where $\mathcal{K}={\left\{\left({a}_k,{b}_k\right)\right\}}_{k=1,\dots,\ K}$ is the support of the bivariate variable (A,B). The joint distribution F _A,B (a,b) is defined by the probabilities p _k  = ℙ(A = a _k , B = b _k  ) _{k = 1,…, K} of all the K admissible couples of values (a _k ,b _k ). 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 $\mathcal{K}$:

$$ rmstD\left({t}^{\ast}\right)={\displaystyle {\sum}_{k=1}^K{p}_k{rmstD}_k\left({t}^{\ast}\right)}, $$

(19)

with the conditional restricted mean survival time ${rmstD}_k\left({t}^{\ast}\right)$ defined for a couple (a _k, b _k ) as:

$$ rmst{D}_k\left({t}^{*}\right)={\displaystyle {\int}_0^{t^{*}} \exp \left(-\frac{ \log (2)}{5}\cdot t\cdot {e}^{a_k+1/2\left(\beta +{b}_k\right)}\ \right)dt-{\displaystyle {\int}_0^{t^{*}} \exp \left(-\frac{ \log (2)}{5}\cdot t\cdot {e}^{a_k-1/2\left(\beta +{b}_k\right)}\ \right)dt}} $$

(20)

Simulation scenarios

In different scenarios we varied the strength of between-trial heterogeneity for the baseline hazard (low with σ ² = 0.01 and high with σ ² = 0.10) and for the treatment effect (τ ² = 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}^{\ast }$= 5 years and ${t}^{\ast }$= 10 years. These two values were chosen to illustrate a scenario in which all trials still have patients at risk at ${t}^{\ast }$ (5 years) and a scenario in which some trials’ follow-up are shorter than the time of restriction (${t}^{\ast }$= 10 years). The average administrative censoring rate ranged across scenarios from 49 to 52 % at ${t}^{\ast }$= 5 years and from 38 to 40 % at ${t}^{\ast }$= 10 years. In the case of no overall treatment effect at all (β = 0, τ ² = 0) and no baseline heterogeneity (σ ² = 0), the restricted mean survival time was equal to 3.6 years at ${t}^{\ast }$= 5 years and 5.4 years at ${t}^{\ast }$= 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\left({t}^{\ast}\right)$ minus the true value; the empirical standard error (ESE), defined as the standard deviation of the $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$, notably on the long term (${t}^{\ast }$= 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).

Table 1 Simulation results for comparisons of methods in estimating the difference in restricted mean survival time. Scenario with 5 trials and 200 patients per trial and with proportional hazards

Fig. 1

Graphical comparison at 5 years (panel a) and at 10 years (panel b) of methods in terms of bias, empirical and average standard error with proportional hazards; 5 trials and 200 patients per trial; β = −0.7. Black squares represent the average bias of the rmstD estimated by each method for a particular scenario. Black and purple vertical lines represent the 95 % confidence intervals of the bias based on respectively the empirical and average standard errors. The horizontal red line indicates the absence of bias in the rmstD estimation. β: Size of treatment effect (=log(HR)); σ²: baseline hazard heterogeneity; τ²: treatment effect heterogeneity; CI: Confidence interval; rmstD: difference in restricted mean survival time

In scenarios with higher treatment effect heterogeneity (τ ² = 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 τ ² = 0.10. When varying the baseline hazard heterogeneity between trials, no relevant impact was noted neither on the bias nor on the standard error.

Fig. 2

Graphical comparison at 5 years (panel a) and at 10 years (panel b) of methods in terms of bias, empirical and average standard error with non-proportional hazards; 5 trials and 200 patients per trial; β = −0.7. Black squares represent the average bias of the rmstD estimated by each method for a particular scenario. Black and purple vertical lines represent the 95 % confidence intervals of the bias based on respectively the empirical and average standard errors. The horizontal red line indicates the absence of bias in the rmstD estimation. β: Size of treatment effect (=log(HR)); σ²: baseline hazard heterogeneity; τ²: treatment effect heterogeneity; CI: Confidence interval; rmstD: difference in restricted mean survival time

Table 2 Simulation results for comparisons of methods in estimating the difference in restricted mean survival time. Scenario with 5 trials and 200 patients per trial, and with non-proportional hazards

With both proportional and non-proportional hazards, for β = −0.7, the Peto-quintile method showed a bias which was negligible at ${t}^{\ast }$= 5 years but much higher at ${t}^{\ast }$= 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\left({t}^{\ast },=,5\right)$ = −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\left({t}^{\ast },=,10\right)$ = 0.30 years (Fig. 2-b).

In terms of standard error, lower values were found for both ESE and ASE at ${t}^{\ast }$ = 5 years (Fig. 1-a and Fig. 2-a) as compared at ${t}^{\ast }$= 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 (σ ² = τ ² = 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 (τ ² = 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\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$, we selected ${t}^{\ast }$= 5 years and ${t}^{\ast }$= 10 years, as these were the two time points of clinical interest in the publications of MAC-NPC and MAC-NPC2.

Meta-Analysis of Chemotherapy in Nasopharynx Carcinoma (MAC-NPC)

The data from the MAC-NPC [34] included 1,975 patients in 11 treatment comparisons. The pooled HR estimated with a fixed effect model was 0.82 (95 % CI: [0.71;0.94]), indicating a significant improvement in overall survival with RT plus CT (p = 0.006). The treatment effect heterogeneity was significant (Q-test: p = 0.03; Higgins’ I² = 50 %) which was explained by the timing of CT. The pooled HR estimated with a DerSimonian-Laird random effects model [28] was 0.82 (95 % CI: [0.66;1.02], p = 0.08).

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\left({t}^{\ast}\right)$ ranged from 0.17 to 0.23 years at ${t}^{\ast }$= 5 years and from 0.46 to 0.55 years at ${t}^{\ast }$= 10 years across the estimation methods (Table 3). For Pooled Kaplan-Meier and Pooled Exponential using a random effects model, the $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$. 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\left({t}^{\ast}\right)$ and SE ($rmstD\left({t}^{\ast}\right)$) were found at ${t}^{\ast }$= 5 years as compared to at ${t}^{\ast }$= 10 years. Figure 3 displays the forest plot for trial-specific ${rmstD}_j\left({t}^{\ast}\right)$ and overall $rmstD\left({t}^{\ast}\right)$ estimated using Pooled Kaplan-Meier with DerSimonian-Laird random effects at ${t}^{\ast }$= 10 for the MAC-NPC meta-analysis. Figure 4 displays the overall rmstD(t ^*) estimated by Pooled Kaplan-Meier with DerSimonian-Laird random effects when varying ${t}^{\ast }$; it shows that the $rmstD\left({t}^{\ast}\right)$ is not significantly different from 0 for ${t}^{\ast }$ϵ [0;10] years. The same graphic for the overall $rmstD\left({t}^{\ast}\right)$ is displayed as Additional file 1: Figure S1.

Table 3 Results for comparisons of methods in estimating the difference in restricted mean survival time (rmstD) in MAC-NPC and MAC-NPC2 meta-analyses

Fig. 3

Forest plot for differences in restricted mean survival time estimated at 10 years using the Pooled Kaplan-Meier method with random effects applied to the MAC-NPC meta-analysis. Each trial is represented by a square, the center of which denotes the difference in restricted mean survival time (rmstD) for that trial comparison, with the horizontal lines showing the 95 % confidence intervals (CI). The size of the square is directly proportional to the amount of information contributed by the trial. The diamond represents the overall rmstD, with the center denoting the rmstD and the extremities the 95 % CI. The rmstDs are expressed in year

Fig. 4

Difference in restricted mean survival time estimated using the Pooled Kaplan-Meier method with random effects as a function of the time horizon t ^* in the MAC-NPC meta-analysis. The solid black line represents the rmstD(t ^*) plotted as a function of the horizon t ^*. The dashed horizontal line indicates the absence of a treatment effect (rmstD(t ^*) = 0). The grey area corresponds to the pointwise 95 % confidence interval. MAC-NPC: Meta-Analysis of Chemotherapy in Nasopharynx Carcinoma; rmstD: difference in restricted mean survival time

Update of Meta-Analysis of Chemotherapy in Nasopharynx Carcinoma (MAC-NPC2)

The MAC-NPC2 [35], the update of the MAC-NPC, included new trials as well as updated follow-up for trials included in the MAC-NPC (N = 5,028 patients within 23 comparisons). For overall survival, a significant pooled HR of 0.79 (95 % CI: [0.73;0.86], p < 0.001) in favor of CT + RT was obtained with a fixed effect model. In the MAC-NPC2, there was less evidence of treatment effect heterogeneity (Q-test: p = 0.09; Higgins’ I² = 30 %) than in the MAC-NPC. The pooled HR with a DerSimonian-Laird random effects model [28] was 0.79 (95 % CI: [0.70;0.87], p <0.001).

The pooled p-value test (p = 0.16) suggested that the overall proportional hazards assumption was appropriate. The $rmstD\left({t}^{\ast}\right)$ ranged from 0.17 to 0.21 years at ${t}^{\ast }$= 5 years and from 0.53 to 0.59 years at ${t}^{\ast }$= 10 across the estimation methods (Table 3). The $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ was lower in the MAC-NPC2 with a $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ and the overall $rmstD\left({t}^{\ast}\right)$ at ${t}^{\ast }$= 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\left({t}^{\ast}\right)$) 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 more intuitive interpretation [3, 5, 6]. Furthermore, the $rmstD\left({t}^{\ast}\right)$ is directly related to cost-effectiveness analysis as it is the denominator of the incremental cost-effectiveness ratio, so one can use the $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ from IPD meta-analysis is available to date [6].

In this study, we compared different methods to estimate the $rmstD\left({t}^{\ast}\right)$ 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}^{\ast }$= 10 years and an even larger bias at t ^* = 5 years. Peto-quintile underestimated the $rmstD\left({t}^{\ast}\right)$, except with non-proportional hazards at ${t}^{\ast }$= 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 random effects underestimated the standard error of the $rmstD\left({t}^{\ast}\right)$. Overall, the Pooled Kaplan-Meier method with DerSimonian-Laird random effects formed the best compromise in terms of bias and variance for estimating the $rmstD\left({t}^{\ast}\right)$ from IPD meta-analysis.

In the IPD meta-analyses studying the effect of chemotherapy (CT) plus radiotherapy (RT) versus RT alone in nasopharynx carcinoma, the $rmstD\left({t}^{\ast },=,10,\kern0.5em ,\mathrm{years}\right)$ estimated using the Pooled Kaplan-Meier method with DerSimonian-Laird random effects was 0.49 years (95 % CI: [−0.06;1.03], p = 0.08) in the MAC-NPC [34] and 0.59 years (95 % CI: [0.34;0.84], p < 0.0001) in its updated version MAC-NPC2 [35]. In other words, the addition of CT to RT extended the 10-year mean survival time by 7.1 months (95 % CI 4.1;10.1) in MAC-NPC2. We believe the clinical interpretation with the $rmstD\left({t}^{\ast}\right)$ is more intuitive than the one derived from the pooled hazard ratio with DerSimonian-Laird random effects of 0.79 (95 % CI 0.70−0.87) in MAC-NPC2.

The $rmstD\left({t}^{\ast}\right)$ 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\left({t}^{\ast}\right)$ 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 meta-analyses in nasopharynx carcinoma, the heterogeneity was slightly higher when pooling the ${rmstD}_j\left({t}^{\ast}\right)$ 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² = 50 % for HR as compared to p = 0.02, I² = 54 % for $rmstD\left({t}^{\ast },=,10,\kern0.5em ,\mathrm{years}\right)$ (Fig. 3). For the MAC-NPC2, this increase was more pronounced with p = 0.09, I² = 30 % for HR versus p = 0.01, I² = 45 % for $rmstD\left({t}^{\ast },=,10,\kern0.5em ,\mathrm{years}\right)$ (Additional file 1: Figure S1). Wei and colleagues showed a similar trend in their second example (p = 0.47, I² = 0 % for HR and p = 0.20, I² = 24 % for $rmstD\left({t}^{\ast },=,5,\kern0.5em ,\mathrm{years}\right)$) [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 ) = σ ² = 0 and Var(B _j ) = τ ² = 0. Testing for the presence of treatment effect heterogeneity (τ ² = 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 (σ ² = 0) in proportional hazard models. Rondeau et al. [19] tested both the baseline hazard (σ ² = 0) and the treatment effect heterogeneity between trials (τ ² = 0) using a mixture of χ ² 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\left({t}^{\ast}\right)$ even for aggregate data. However, we suggest that clinical publications for single (multicenter) clinical trial or IPD meta-analysis should report the $rmstD\left({t}^{\ast}\right)$ at different time horizons ${t}^{\ast }$ of clinical interest in addition to the hazard ratio. This way the $rmstD\left({t}^{\ast}\right)$ 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 non-parametric 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\left({t}^{\ast}\right)$. In addition, in this simulation study, we only considered balanced trials and we did not vary the administrative censoring rate.

The $rmstD\left({t}^{\ast}\right)$ is inherently dependent on the choice of ${t}^{\ast }$. Also, we showed that its standard error gets larger as ${t}^{\ast }$increases (Fig. 1). Karrison recommended to choose a maximum time horizon ${t}^{\ast }$ such that SE(S(${t}^{\ast }$)) 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}^{\ast }$ that censor a lot of data with a high loss of information, and high values of ${t}^{\ast }$ that need a massive use of extrapolation. Wei and colleagues stated that the choice of ${t}^{\ast }$ should also be of clinical interest, and they suggested plotting the $rmstD\left({t}^{\ast}\right)$ against ${t}^{\ast }$ to see how the treatment effect varies over time. In MAC-NPC for instance such a plot shows that the $rmstD\left({t}^{\ast}\right)$ was not significantly different from 0 with ${t}^{\ast }$ϵ [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\left({t}^{\ast}\right)$ 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}^{\ast }$, some trials included in the IPD meta-analysis may have a follow-up not long enough to reach ${t}^{\ast }$. In our study, for such trials, we used the extrapolation method proposed by Brown et al. [32] until ${t}^{\ast }$ 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\left({t}^{\ast}\right)$) 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\left({t}^{\ast}\right)$ 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.