Bayesian survival extrapolation for cost-effectiveness analysis: a case study of RELAY for ramucirumab in combination with erlotinib in the treatment of non-small-cell lung cancer

Abstract

Aim

Increasing trend for progression-free survival (PFS)-based primary endpoint in oncology has led to lack of mature overall survival (OS) data at the time of approval. To address this evidence gap in economic evaluations, we used a joint Bayesian approach to predict survival outcomes using immature OS data from the RELAY trial.

Methods

Patient data from RELAY and systematic literature review (SLR) of phase 3 randomized clinical trials with hazard ratio (HR) estimates of mature PFS and immature OS were considered. OS and PFS were analyzed individually using a univariate model; bivariate analysis was performed using a joint model based on modified Bayesian normal induced copula estimation model. First, a Bayesian univariate model incorporated informative priors based on predicted HR and acceleration factor for OS and PFS. Second, a Bayesian-based joint model of RELAY PFS and OS data was based on the correlation between PFS and OS established in trials of similar populations. Marginal distribution of PFS was used to estimate the same for OS.

Results

Publications (N = 122) of first-line treatments in patients with epidermal growth factor receptor (EGFR)-mutated non-small cell lung cancer were identified in the SLR, of which 36 trials were linked to RELAY. Twenty-six trials with HR data were used. The univariate model could predict OS with reduced uncertainty compared with the frequentist approach. In the joint model, the marginal OS distribution borrowed strength from the marginal PFS distribution through the established correlation coefficient.

Limitations

Bayesian approach was successfully used in RELAY analysis but may not be universally applied to oncology trials due to the different associations of OS and PFS and different trial patient populations.

Conclusions

We demonstrated that both the univariate and joint Bayesian models reduced uncertainty in predicting OS compared to frequentist method. The methodology introduced here will have potential applications in clinical decision-making for other oncology trials.

Keywords:

• Overall survival
• progression-free survival
• health technology assessment bodies
• Bayesian modeling
• RELAY
• ramucirumab

JEL CLASSIFICATION CODES:

• C11
• C53
• C44
• C6

Introduction

In the field of oncology, the criteria by which drugs receive regulatory approval are ever evolving. An increasing number of drug approvals in the metastatic setting are now based primarily on progression-free survival (PFS) outcomes¹. Correspondingly, the increase in PFS-based drug approvals, along with an increase in expedited approvals of anticancer and biologics in drug development, has resulted in the lack of mature overall survival (OS) data at the time of regulatory approval.

The immaturity of OS data presents a significant challenge for health technology assessment (HTA) bodies in countries where a cost-effectiveness model is required for reimbursement purposes. The OS data were used to calculate life years, and quality adjusted life years in cost-effectiveness analyses. Therefore, mature OS data are required to analyze drug cost-effectiveness.

In the countries where cost-effectiveness analysis is critical for decision making, uncertainty in OS outcomes is a significant challenge when evaluating a new therapy. Partitioned survival models in a frequentist framework have become the standard in recent oncology HTA submissions in the metastatic setting. However, these model results are highly sensitive to immature OS data. In cases where patients are living longer, OS data are often immature, leading to uncertainty associated with survival outcomes and reluctancy of HTA bodies to recommend². To this end, due to the lack of mature OS evidence, many promising anticancer therapies with superior PFS results may not become fully accessible to patients soon after approval. Thus, there is a pressing need to address payer agencies’ criticism about uncertainty resulting from immature OS in economic models.

Non-small cell lung cancer (NSCLC) represents around 85% of all lung cancers^3–5. Each year, approximately 2.1 million individuals are diagnosed with lung cancer; and the majority of those with NSCLC present with advanced or metastatic disease states^6–9. The finding that, genetically, NSCLC is a highly heterogenous disease, and that these genetic alterations that drive initial tumor growth and progression in some forms of NSCLC can be targeted, has redefined the treatment landscape of lung cancer¹⁰. In the United States, it is estimated that around 19% of people diagnosed with NSCLC have an epidermal growth factor receptor (EGFR) mutation¹¹. For these patients, EGFR tyrosine kinase inhibitors (TKIs) are the first-line standard of care¹². However, although patients respond positively to TKIs at first, resistance is inevitable and usually develops it in under 19 months of treatment^13,14. The use of combination therapy regimens to further slow NSCLC disease progression and resistance to TKIs is a growing area of interest. The combination of vascular endothelial growth factor (VEGF) inhibitors with EGFR TKIs is one such approach, which is already proving promising.

The randomized phase 3 RELAY trial (NCT02411448) was conducted to investigate the efficacy of ramucirumab, a selective antagonist of VEGF receptor 2 (VEGFR2), in combination with erlotinib versus placebo and erlotinib, as first-line treatment in adults with metastatic EGFR-mutated NSCLC. At the time of primary analysis, median PFS was significantly longer in the ramucirumab plus erlotinib group (19.4 months (95% CI 15.4–21.6)) compared to the placebo plus erlotinib group (12.4 months (95% CI 11.0–13.5)), with a stratified hazard ratio (HR) of 0.59 (95% CI 0.46–0.76; p < .0001)¹⁵. Based on a statistically significant and clinically meaningful improvement in PFS, ramucirumab plus erlotinib received regulatory approval^16,17.

In countries where reimbursement decisions are based on HTA bodies, the cost-effectiveness of ramucirumab plus erlotinib will be compared to other treatment options. At the time of the primary analysis of RELAY, interim OS demonstrated no detriment¹⁵. However, HTA bodies have previously criticized the immaturity of OS data for first-line treatments in patients with EGFR-mutated NSCLC, including osimertinib (FLAURA trial − 25% mature) and dacomitinib (ARCHER trial – >45% mature), despite PFS data being fully mature, statistically significant, and clinically meaningful at the time of the primary analysis^18,19. In the ramucirumab RELAY trial, the OS data maturity level was 26% (79 of 300 planned OS events at final analysis with censoring rate 82% (370 of 449))¹⁵.

Using the ramucirumab RELAY trial as a case study, our modeling framework addresses the issue of immature OS data by using two Bayesian approaches to improve projected OS estimates for patients with NSCLC within the study period and the extrapolated survival probability beyond the study period from immature data and external data. First, a univariate modeling approach was applied, incorporating informative priors to inform RELAY OS marginal distribution. External data, including HR, Kaplan–Meier (KM), and accelerated failure time (AFT) coefficient estimates for OS and PFS derived from 26 other trials in patients with first-line treatment naïve patients with EGFR-positive NSCLC, were used to formulate the informative prior. By combining these pre-existing data, we reveal that treatment effect in PFS can inform treatment effect in OS. Due to this correlation, a predicted effect in OS (HR and AFT) was derived for RELAY using the observed effect in PFS of the trial, and the resulting predicted effect in OS (HR and AFT) was used as the informative prior in a Bayesian survival model.

The second approach is a Bayesian-based joint model of RELAY PFS and OS data. Depending on the extent to which RELAY PFS is associated with RELAY OS, the marginal OS distribution is expected to borrow strength from the marginal distribution of PFS through the correlation allowed by the joint model. In addition, the correlation coefficient between OS and PFS will be informed by a prior distribution constructed by data from trials conducted in a similar population to the RELAY trial. Here, we applied this joint Bayesian approach to predict RELAY survival outcomes based on the immature OS data.

Methods

Sources of data

Two main sources of data were used in this analysis.

Individual patient data

Individual patient data (IPD) from the primary analysis of the randomized, placebo-controlled, phase 3 RELAY trial (NCT02411448) are one of the main data sources (data cutoff occurred on 23 January 2019). This included mature PFS, however, OS was immature (censoring rate of 82%) at the time of the primary analysis. Final OS, a secondary objective of RELAY and that the study is not powered to detect an OS difference, will be evaluated after approximately 300 deaths have occurred. RELAY investigated the efficacy of ramucirumab in combination with erlotinib in first-line treatment of adults with metastatic NSCLC, versus placebo plus erlotinib.

Systematic literature review data

The other source of data was a systematic literature review (SLR) conducted via Ovid. A SLR was performed to identify randomized phase 3 trials in first-line patients with advanced (stage III/IV) EGFR-mutated NSCLC, to the RELAY trial. Studies with available HR estimates of PFS and OS were used for this analysis²⁰.

Analysis data type

OS analysis – univariate model

First, a linear predictive model was developed using data from 26 randomized phase 3 trials of first-line treatments for patients with EGFR-mutated NSCLC with OS HR (on log scale) acting as a dependent variable and PFS HR (on log scale) as a predictor. This predicted OS HR (log scale) for RELAY was incorporated in the Bayesian model as the mean of the prior distribution of log HR assumed to be normally distributed for parametric distributions allowing proportional hazards (exponential, Weibull and Gompertz). $$\begin{array}{l}\log \left({\text{HR}}_{\text{os}}\right)=\alpha +\beta .log\left({\text{HR}}_{\text{pfs}}\right).\\ \text{The predicted OS HR using the RELAY PFS was denoted by }\\ {\widehat{HR}}_{\text{OS}}^{\text{RELAY}}.\end{array}$$

A similar linear model based on treatment effect derived from AFT was considered for the distributions that do not have proportional hazards (log-logistic, log-normal, and generalized gamma). KM curves were published for 17 trials only. These curves were extracted and pseudo IPD were reconstructed based on the method described by Guyot et al²¹. The resulting survival data were used to construct the linear regression model that predicts the OS AFT based on PFS AFT. Subsequently, the OS acceleration factor for RELAY was predicted using the observed RELAY acceleration factor of PFS for each of the AFT distributions.

PFS analysis – univariate model

A univariate Bayesian model assuming non-informative prior distributions was considered for PFS. The posterior distribution was based solely on the observed RELAY PFS data. The resulting univariate model PFS should be comparable to the RELAY trial PFS results, which were analyzed using a frequentist approach.

OS and PFS analysis – joint model

The joint model is based on a Bayesian normal induced copula estimation model modified to analyze five additional distributions beyond the exponential distribution that was reported by Fu et al²². This model assumed a known correlation coefficient value between OS and time-to-progression (TTP). As TTP data were not routinely reported in the identified trials from the SLR, the correlation coefficient between OS and PFS was considered. The estimation of the correlation coefficient also required an assumption about the correlation between OS and PFS. Two approaches for estimating the Pearson correlation were considered using 17 randomized trials of first-line treatments for patients with EGFR-mutated NSCLC, where the HR estimates were available for both OS and PFS. The first approach was based on Fleischer et al²³. assuming exponential distribution for OS and PFS. The second was based on the finding of Li and Zhang²⁴ with the assumption that OS and PFS follow the Weibull distribution.

Various distributions were considered for the likelihood function, including the six distributions of interest.

Model for treatment effect

Treatment effect is modeled as the acceleration factor in the accelerated failure models (log-logistic, log-normal, and generalized gamma) and as a constant HR for distributions that allow proportional hazards (exponential, Weibull and Gompertz).

Bayesian modeling approach

OS analysis – univariate model

For the univariate Bayesian modeling approach, non-informative prior distributions were used for all parameters of the six parametric distributions except for treatment effect. Non-informative prior distributions may be eventually extended to minimal information in the priors. An informative prior was used for the treatment effect measured by log HR or log AFT. Specifically, the coefficient followed a normal distribution with mean and variance corresponding to the predicted OS effect for RELAY ${\widehat{HR}}_{\text{OS}}^{\text{RELAY}}$ (transformed on the log scale) and its variance, respectively. The formulation was given by the following equation: $$\log \left({\text{HR}}_{\text{os}}\right)\backsim \mathit{normal} {(\log (\widehat{HR}}_{\text{OS}}^{\text{RELAY}}), \mathit{var}({\log (\widehat{HR}}_{\text{OS}}^{\text{RELAY}}\left)\right)$$

A similar linear model based on treatment effect derived from AFT was considered for the distributions that do not rely on proportional hazards and HRs (log-logistic, log-normal, and generalized gamma). Only the 17 trials that report that KM were used to predict the OS acceleration factor for the RELAY using the acceleration factor of PFS. $$\log \left({\text{AFT}}_{\text{os}}\right)\backsim \mathit{normal}\mathrm{ }{(\log (\widehat{\mathit{AFT}}}_{\text{OS}}^{\text{RELAY}}),\mathrm{ }\mathit{var}({\log (\widehat{\mathit{AFT}}}_{\text{OS}}^{\text{RELAY}}\left)\right)$$

PFS analysis – univariate model

Non-informative prior distributions were used for parameters of the six parametric distributions.

OS and PFS analysis – joint model

A uniformly distributed prior for correlation coefficient was considered for the joint modeling Bayesian approach, the lower and upper bound parameters correspond respectively to minimum and maximum estimated correlations among the 17 studies in the SLR from which the estimation was possible: $$\rho \backsim \mathit{uniform}\left(\min \left({\widehat{\rho }}^k\right),\max \left({\widehat{\rho }}^k\right)\right) \mathrm{for}\mathrm{ }k=1,2,\dots ,17.$$

${\stackrel{\mathrm{`}}{\rho }}^k$ is the estimated Pearson correlation for the trial k from the SLR results.

Likelihood function

OS analysis – univariate model

Bayesian methods involved formal combination of a prior probability distribution with the likelihood obtained from a sampling distribution used to model the observed data to obtain a posterior probability distribution of model parameters. The likelihood informed us about the extent to which different values for the parameter of interest were supported by the data.

The RELAY studied IPD were right censored for OS. The $f_{\theta }\left(t\right)$ and $S_{\theta }\left(t\right)$ were the probability density function and the survival function of the parametric distribution, respectively, with $\theta$ as parameter vector. Each patient then contributed $f_{\theta }\left(t_j\right)$ at observed time $t_j$ when the event was death and contributed $S_{\theta }\left(t_j\right)$ when the event was censoring.

The probability density function, $f_{\theta }\mathrm{(}\text{.}\mathrm{)}$ and the survival function $S_{\theta }\mathrm{(}\text{.}\mathrm{)}$ were known for all six distributions of interest: exponential, Weibull, Gompertz, log-logistic, log-normal, and generalized gamma. The data consist of $n$ replications denoted as $\mathrm{(}{X}_i,{\Delta }_i\mathrm{)}$ for patient $i=1,\dots ,n,$ where $X_i$ was the follow-up time for OS of patient $i,$ which was time to death or censoring whichever happened first and ${\Delta }_i$ was binary variable equal to 1 when the patient $i$ died and 0 when the patient was censored. Assuming non-informative censoring, the likelihood function was: $$L\left(\theta |\text{Data}\right)={\prod }_{i=1}^n{f_{\theta }\left(t_i\right)}^{{\Delta }_i}{{\mathrm{S}}_{\theta }\left(t\right)}^{(1- {\Delta }_i)}$$

PFS analysis – univariate model

The univariate model for PFS analysis used the same likelihood function formulation as described in the univariate model analysis of OS.

OS and PFS analysis – joint model

The likelihood function was specified in the same way as the method proposed by Fu et al²². The data consisted of $n$ replications denoted as $\mathrm{(}{X}_i,{\mathrm{\Delta }}_i^1,Y_i,{\mathrm{\Delta }}_i^2\mathrm{)}\text{, }i=1,\dots ,n,$ where $X_i$ was the follow-up time for OS which was time to death or censoring whichever happened first and ${\mathrm{\Delta }}_i^1$ equals 1 when the patient died and 0 when the patient was censored. Similarly, $Y_i$ was the follow-up time for PFS and ${\mathrm{\Delta }}_i^2$ was the corresponding status.

There are four mutually exclusive cases for each patient:

1. Case 1: OS (death) was observed (${\mathrm{\Delta }}_i^1=1$) and TTP (progression) was also observed (${\mathrm{\Delta }}_i^2=1$)

The contribution to the likelihood function was in this case determined by the probability density function of the joint distribution (OS, TTP).

• 2. Case 2: OS (death) was observed (${\mathrm{\Delta }}_i^1=1$) but TTP (progression) was not observed (${\mathrm{\Delta }}_i^2=0$)

The contribution to the likelihood function was in this case determined by the joint probability that death was observed and TTP was censored by death.

• 3. Case 3: OS time was censored at time ($t_{i1},$ ${\mathrm{\Delta }}_i^1=0$) though TTP (progression) was observed (${\mathrm{\Delta }}_i^2=1$)

The contribution to the likelihood function was in this case determined by the joint probability that TTP was observed and death was censored.

• 4. Case 4: OS and TTP times were censored at time $t_i,$ in other words ${\mathrm{\Delta }}_i^1=0$ and ${\mathrm{\Delta }}_i^2=0.$ The contribution to likelihood function was in this case determined by the joint survival function of OS and TTP (both OS and TTP were censored).

When ${\theta }_1$ and ${\theta }_2$ were the two vectors of parameters in the hazard functions of TTP and OS time, respectively, and $\rho$ was the correlation coefficient between TTP and OS, with non-informative censoring assumption, the likelihood function was given by: $$L\left({\theta }_1,{\theta }_2,\rho |\text{Data}\right)={\prod }_{i=1}^n{f\left(t_{i1},t_{i2}\right)}^{{\mathrm{\Delta }}_i^1{\mathrm{\Delta }}_i^2}{S\left(t_{i1},t_{i2}\right)}^{{\mathrm{(}\text{1}\text{-}\mathrm{\Delta }}_i^1\left){\mathrm{(}\text{1}\text{-}\mathrm{\Delta }}_i^2\right)}{f_1\left(t_{i1},t_{i2}\right)}^{{\mathrm{\Delta }}_i^1{\mathrm{(}\text{1}\text{-}\mathrm{\Delta }}_i^2)}{f_2\left(t_{i1},t_{i2}\right)}^{{\mathrm{(}\text{1}\text{-}\mathrm{\Delta }}_i^1){\mathrm{\Delta }}_i^2}$$ where $f$ is the probability density function for joint distribution of TTP and OS; and $f_1$ and $f_2$ were the partial derivatives of joint distributions relative to TTP and OS, respectively.

A similar approach to Fu et al²². based on copula was taken to generate the joint distribution of OS and TTP using the observed Pearson correlation between OS and TTP from the SLR studies. In addition to exponential distribution used in Fu et al²²., the five parametric distributions usually required by HTA bodies in frequentist context were fitted to. Thus, the following six distributions were shown in the Bayesian joint model of OS and TTP: exponential, Weibull, Gompertz, log-logistic, log-normal, and generalized gamma. Afterwards, the PFS was derived based on OS and TTP by simulation as proposed by Fu et al²².

Results

Results of the systematic literature review

In total, 122 publications of first-line treatments in patients with EGFR-mutated NSCLC were identified in SLR, of which 57 were found to be unique trials investigated in a first-line setting²⁰. Thirty-six randomized controlled trials reported data for any endpoint assessed and directly or indirectly connected to a treatment in the RELAY trial through TKI monotherapy or combination therapies, as previously described in the meta-analysis by Haeussler et al²⁰. However, for this analysis, trials identified during the literature review that did not report HR for both PFS and OS were excluded as both data points were required to develop the linear model. Based on these criteria, 26 of the identified 36 trials were included in this study. The baseline characteristics of the target population varied across the included trials. While the distribution of age, sex, percentage of ever smokers, and ECOG PS were widely varied, the distribution of NSCLC histological subtype was consistent across the trials, with the exception of TORCH trial²⁰. Overall, the similarity assumption held except for patients with brain and central nervous system (CNS) metastases. As previously demonstrated, through the application of unrelated means models to the data, the assumptions of consistency and transitivity were deemed to hold. This indicated that there were comparable distributions of known effect modifiers between studies²⁰.

OS analysis – univariate model

The interim stratified OS HR 0.83 (95% CI 0.53–1.30) was evaluated at the time of data lock¹⁵ out of the 26% of planned events for final OS analysis. First, a univariate Bayesian modeling approach was applied to extrapolate OS using informative prior to inform the marginal distribution from 26 trials similar to the RELAY trial. The fitted Weibull, exponential, and Gompertz curves are plotted against the nonparametric (KM) curves in Figure 1, Supplementary Figures 1 and 2, respectively, from which an adequacy of fit and an increasing separation between the curves over time was observed. The resulting estimates of model parameters were consistent with the frequentist immature OS; however, the credible interval limits were much smaller than the 95% CI associated with the frequentist estimates. Upon the final analysis at 300 total OS events, the projected KM estimate of median OS was predicted to be longer for the ramucirumab plus erlotinib arm compared to the placebo plus erlotinib arm (Table 1). The resulting HR for OS obtained with the Bayesian univariate model with informative prior was 0.82 (95% credible interval 0.75–0.90), which is consistent with the point estimate of the frequentist HR; however, the 95% credible interval was smaller than the 95% CI associated with the frequentist HR (Table 2).

Figure 1. Placebo plus erlotinib and ramucirumab plus erlotinib Kaplan–Meier’s survival curve with estimated survival probability via the Weibull distribution with informative prior. Abbreviations. Erl, Erlotinib; PBO, Placebo; OS, Overall survival; RAM, Ramucirumab.

Table 1. Summary of posterior distribution of parameters for overall survival (OS) using informative prior distribution.

	Exponential
	Intercept	Treatment effect
Mean	4.694	0.194
Median (95% CrI)	4.690 (4.478–4.930)	0.194 (0.105–0.282)
	Weibull
	Intercept	Treatment effect	Scale
Mean	4.130	0.185	0.652
Median (95% CrI)	4.121 (3.903–4.41)	0.184 (0.097–0.272)	0.648 (0.534–0.795)
	Gompertz
	Intercept	Treatment effect	Gamma
Mean	5.394	0.194	0.051
Median (95% CrI)	5.385 (4.96–5.869)	0.193 (0.107–0.284)	0.051 (0.025–0.078)
	Log-logistic
	Intercept	Treatment effect	Gamma
Mean	6.736	0.169	1.680
Median (95% CrI)	6.709 (5.759–7.855)	0.169 (0.06–0.279)	1.672 (1.373–2.035)
	Log-normal
	Intercept	Treatment effect	Gamma
Mean	4.330	0.002	1.283
Median (95% CrI)	4.319 (4.015–4.71)	0.002 (–0.345 to 0.344)	1.274 (1.074–1.539)
	Generalized gamma
	Intercept	Treatment effect	Scale	Shape
Mean	18.324	−1.713	18.788	−4.802
Median (95% CrI)	18.247 (6.932–31.146)	−1.751 (–5.348 to 2.205)	19.209 (10.029–25.414)	−4.362 (–10.503 to 1.862)

Table 2. Hazard ratio for overall survival using a Bayesian and frequentist approach.

	Bayesian
RAM + ERL vs. PBO + ERL	HR (95% credible interval)
	0.82 (0.75, 0.90)
	Frequentist
RAM + ERL vs. PBO + ERL	HR (95% Wald confidence limit)
	0.83 (0.53, 1.30)

PFS analysis – univariate model

Using a frequentist approach to predict PFS, it was found that ramucirumab plus erlotinib group (median PFS, 19.4 months (95% CI 15.4–21.6)) PFS was significantly longer than placebo plus erlotinib group (median PFS, 12.4 months (95% CI 11.0–13.5)), with a stratified HR of 0.59 (95% CI 0.46–0.76)¹⁵. RELAY PFS data were re-analyzed using a Bayesian framework to ensure consistency between the frequentist and Bayesian approaches (Supplementary Table 1). In this case, a univariate Bayesian model with a non-informative prior distribution was considered for PFS; the posterior distribution was based solely on the observed data in the RELAY trial. The resulting point and interval estimates were expected to be the same for both approaches; the Bayesian framework provided additional information about model parameter distributions. Near equivalent log HRs for PFS were obtained with the Bayesian univariate model with vague Gaussian non-informative prior (HR 0.59; 95% CI 0.46–0.75), and the parametric frequentist model (HR 0.59; 95% CI 0.46–0.76) (Figure 2 and Supplementary Table 2). The fitted exponential and Gompertz curves were plotted over the PFS KM curves in Supplementary Figures 3 and 4, respectively.

Figure 2. Placebo plus erlotinib and ramucirumab plus erlotinib Kaplan–Meier’s survival curve with estimated progression-free survival probability via the Weibull distribution with informative prior. Abbreviations. Erl, Erlotinib; PBO, Placebo; PFS, Progression-free survival; RAM, Ramucirumab.

After demonstrating that a Bayesian approach could be used to accurately predict PFS using mature data, we next explored whether an OS could be predicted using mature PFS data with less uncertainty.

Model selection

The systematic model selection usually provided in frequentist model fitting was still applicable; the traditional frequentist good-fit statistics akaike information criterion (AIC) and Bayesian information criterion (BIC) for model selection were replaced by the Bayesian model fit statistics deviance information criterion (DIC), Dbar and Dmean and pD. Supplementary Tables 3 and 4 display these statistics for RELAY PFS Bayesian model fitting with non-informative prior.

OS and PFS analysis – joint model

Next, a joint Bayesian modeling approach was applied that assumed a correlation exists between RELAY OS and PFS in patients with NSCLC. The correlation coefficient between OS and PFS was informed by a prior distribution constructed using data from 17 of the 26 randomized trials of first-line treatments in patients with EGFR-mutated NSCLC (where reconstructed KM estimates for both OS and PFS were available). Figure 3 depicts the observed PFS and OS KM curves by treatment arm, overlaid by the estimated Weibull marginal distributions of OS and PFS based on the joint Bayesian model of OS and TTP. The figure demonstrated that the predicted marginal Weibull model for PFS visually fitted well to the observed PFS KM curve. Moreover, the model-estimated curve for OS was aligned with the observed OS data. Comparable results were observed for each distribution (Supplementary Figures 5–8); however, the Weibull distribution was being utilized to illustrate the findings. Table 3 demonstrates that the correlation between OS and PFS (median and mean of posterior) was strong. Interestingly, the correlation coefficients were consistent across each probability distribution.

Figure 3. Fitted Weibull distribution overlaid by Kaplan–Meier’s curves of (A) overall survival, and (B) progression-free survival, by treatment.

Table 3. Summary of posterior distribution of parameters for overall survival (OS) and progression-free survival (PFS) based on the joint model.

	Exponential
	RAM + ERL			PBO + ERL
	OS rate	PFS rate	Correlation	OS rate	PFS rate	Correlation
Mean	0.012	0.052	0.808	0.013	0.072	0.796
SD	0.016	0.005	0.045	0.002	0.006	0.044
Median (95% CrI)	0.012 (0.009–0.015)	0.052 (0.043–0.062)	0.814 (0.703–0.877)	0.0128 (0.010–0.016)	0.072 (0.062–0.084)	0.8021 (0.692–0.866)
	Weibull
	RAM + ERL					PBO + ERL
	OS shape	OS scale	PFS shape	PFS scale	Correlation	OS shape	OS scale	PFS shape	PFS scale	Correlation
Mean	1.460	73.577	1.580	20.334	0.677	1.773	55.871	1.458	15.561	0.642
SD	0.206	15.548	0.117	1.163	0.066	0.233	8.035	0.089	0.865	0.071
Median (95% CrI)	1.4514 (1.085–1.890)	70.772 (51.728–112.115)	1.578 (1.355–1.815)	20.289 (18.157–22.754)	0.682 (0.530–0.790)	1.764 (1.349–2.252)	54.671 (43.813–74.756)	1.456 (1.290–1.638)	15.553 (13.897–17.306)	0.649 (0.490–0.764)
	Gompertz
	RAM + ERL					PBO + ERL
	OS shape	OS rate	PFS shape	PFS rate	Correlation	OS shape	OS rate	PFS shape	PFS rate	Correlation
Mean	0.0376	0.0059	0.0596	0.025	0.718	0.045	0.006	0.036	0.049	0.713
SD	0.0177	0.0017	0.0108	0.005	0.064	0.017	0.002	0.009	0.007	0.065
Median (95% CrI)	0.037 (0.006–0.074)	0.006 (0.003–0.010)	0.060 (0.039–0.081)	0.025 (0.017–0.035)	0.725 (0.578–0.822)	0.045 (0.012–0.079)	0.006 (0.003–0.010)	0.036 (0.017–0.054)	0.049 (0.037–0.063)	0.721 (0.568–0.820)
	Log-Logistic
	RAM + ERL					PBO + ERL
	OS shape	OS scale	PFS shape	PFS scale	Correlation	OS shape	OS scale	PFS shape	PFS scale	Correlation
Mean	1.573	64.382	2.074	15.543	0.668	1.993	48.066	2.111	11.239	0.614
SD	0.218	13.287	0.149	1.039	0.066	0.261	6.743	0.126	0.679	0.070
Median (95% CrI)	1.565 (1.175–2.024)	62.081 (45.672–97.730)	2.071 (1.789–2.374)	15.502 (13.616–17.701)	0.674 (0.526–0.781)	1.982 (1.515–2.529)	47.082 (37.910–63.796)	2.108 (1.870–2.366)	11.222 (9.958–12.627)	0.620 (0.466–0.733)
	Log-normal
	RAM + ERL					PBO + ERL
	OS meanlog	OS sdlog	PFS meanlog	PFS sdlog	Correlation	OS meanlog	OS sdlog	PFS meanlog	PFS sdlog	Correlation
Mean	4.373	1.374	2.687	0.904	0.675	3.988	1.053	2.367	0.868	0.637
SD	0.234	0.173	0.076	0.058	0.065	0.140	0.121	0.068	0.053	0.062
Median (95% CrI)	4.352 (3.978–4.899)	1.359 (1.084–1.763	2.686 (2.543–2.840)	0.900 (0.799–1.025)	0.681 (0.533–0.786)	3.980 (3.763–4.285)	1.026 (0.838–1.317)	2.359 (2.247–2.494)	0.867 (0.786–1.000)	0.643 (0.511–0.756)

Discussion

As a result of the rapid evolution of the oncology drug therapy landscape, in which patients can receive multiple lines of treatment, the design of clinical oncology trials has shifted away from OS toward PFS as the primary endpoint. Thus, OS data are commonly immature during regulatory approval of oncology therapeutics. Consequently, the level of maturity of OS data required for regulatory approval differs from that required for drug appraisal by HTA bodies, creating a discrepancy that needs to be addressed.

While clinical trial data drive regulatory approval, multiple other sources of data can be used in HTA dossiers to inform how interim OS data are likely to mature based on observed PFS. This includes comparative analyses of published PFS and OS data from competitive treatments to inform the relative treatment effects among comparators in cost-effectiveness models. Still, uncertainty around predicting outcomes based on immature OS data in cost-effectiveness models remains a challenge that needs to be addressed to facilitate decision-making by HTA bodies. In this analysis, the Bayesian framework was applied to unify aggregated published data, and patient-level mature PFS and interim OS data from clinical trials in order to reduce uncertainty around outcomes based on immature OS in the cost-effectiveness models.

A standard partitioned survival model in the metastatic setting requires both PFS and mature OS to estimate life years in progression-free and post-progression health states. Furthermore, for surrogate endpoints like PFS to be accepted along with immature OS by HTA bodies, two independent conditions must be met. First, a clear relationship of PFS with OS must be demonstrated, and second, the relationship must then be validated in clinical studies in the relevant disease setting and target population.

Here, our work applied a new unified modeling approach. We used the RELAY trial as a case study to explore alternative solutions for conducting survival analyses in economic modeling with immature OS data. Weibull distribution was utilized to model and depict survival data as both AFT and proportional hazard models can be formulated based on the same. Moreover, Weibull was the best fit model based on deviance information criteria. Equivalent PFS values were obtained using a univariate Bayesian approach without informative prior compared to the traditional frequentist approach.

After determining that a correlation existed in treatment effects between PFS and OS, a prior belief about OS effect could be built using known PFS effect and the established relationship in treatment effects between OS and PFS in other studies of patients with EGFR-mutated NSCLC. It should be noted that the correlation of PFS and OS in RELAY and in other clinical NSCLC studies enables the success of this approach and is an important initial consideration. Our findings demonstrated that the Bayesian univariate model approach could be utilized to predict OS with reduced uncertainty compared to the frequentist approach by borrowing prior information from SLR findings. These findings were similar between the frequentist and univariate Bayesian approaches when no external information was used. For the joint Bayesian model approach, mature PFS data from the RELAY trial were used to inform immature OS data from the trial. The Pearson coefficient established from 17 clinical trials demonstrated that the correlation between PFS and OS was high enough to be used. The correlation coefficient then acted as an informative prior to accurately estimate the marginal distribution of OS. This suggests that, as expected, the marginal OS distribution borrowed strength from the marginal distribution of PFS through the established correlation coefficient. Moreover, the predicted estimate of OS was superior for the ramucirumab plus erlotinib arm compared to the placebo plus erlotinib arm, confirming the interim OS result. While the predictive model does not resolve problems associated with immature OS data, the consistent narrow credible intervals across all distributions reduces the uncertainty.

Based on these findings, the joint Bayesian model is reasonably predictive for OS, and may be a useful tool for economic modeling. As with the traditional frequentist approach, both the univariate and joint Bayesian models allowed fitting of six standard parametric distributions and used the same model selection and validation process per the NICE guidelines. In addition, the Bayesian framework offers the entire posterior distribution for the model parameters that can be directly used in the probabilistic sensitivity analysis (PSA) required by HTA bodies, as an alternative to extensive simulations of clinical trial data within the cost-effectiveness model. The resulting PSA carrying the uncertainty in survival predictions through posterior distributions of parameters may support addressing HTA driven payer agencies’ criticism of immature OS in economic models.

The use of this methodology has the potential to have a positive impact on a variety of stakeholders. From the perspective of HTA bodies, they have the ability to make more informed evaluations of the drug’s cost-effectiveness in cancers where PFS and OS are correlated. Finally, because this approach provides reasonable predictions before mature OS data are seen, medicines may be evaluated and become accessible to patients earlier. In addition, the approach may allow for resources to be further directed at additional clinical development.

Notwithstanding, this approach has limitations in that it cannot be universally applied to all clinical trials because the associations between OS and PFS can differ between cancer types, and it does not account for differences in trial patient populations. In other studies, PFS has been shown to be a poor predictor of OS²⁵. In colorectal cancer, the correlation seems to be well defined, and several authors have shown that PFS is a good predictor of OS. However, the initial approval of bevacizumab for treatment of metastatic breast cancer based on significant PFS results was followed by the disapproval based on insignificant OS results²⁵. Furthermore, some clinical trials may have a limited amount of external data to add to the model’s strength. Thus, the results of the model are a prediction and there is a risk that they may not be the final observed result.

Ultimately, while the proposed methods help reduce the uncertainty around OS in the context of clinical trials, a combination of clinical trial data and real-world data may be ideal for informing this approach, there is a fundamental limitation for all trial data as primary sources for HTA submissions. Both observed KM estimates of survival (OS and PFS) and extrapolations beyond study period would require external validation based real-world data. The proposed approaches in this manuscript are not exceptions in that respect. Nevertheless, we believe that the eligible trials for the construction of the prior correlation must be clinically and methodologically similar to the RELAY trial. In addition, such real-world data are not available for experimental new treatments at the time when the proposed methodology is intended for use. The findings of the proposed methods need to be re-assessed once long-term follow-up data in real-world settings are available. Nevertheless, progression is not assessed systematically in the real-world setting as it is in a clinical trial, and this discrepancy could further introduce potential bias.

Conclusions

As breakthrough cancer therapies continue to extend OS of patients, and especially in the dawning era of precision medicine, we anticipate an increasing number of drug approvals based on strong treatment effects on outcomes such as PFS, accompanied by immature OS data. Here, we demonstrated that both the univariate and joint Bayesian models reduced uncertainty in predicting OS. While Weibull was used to depict the findings, similar results were observed across all distributions analyzed. This work will have a broader impact, with potential applications for other oncology trials in clinical decision-making and optimal treatment decision making when known correlation between PFS and OS exists. If survival prospects are predicted to be positive this may then lead to more informative decision-making among HTA bodies that assess cost-effectiveness of treatments to allow for patient access. In a much wider context, this approach underlines the importance of using a unified framework to bring all available evidence together, including both within and outside the trial, to reduce extrapolation uncertainty of predicting outcomes of clinical trial endpoints. Nevertheless, a combination of clinical trial data and real-world data is ideal for informing this approach for both internal and external validation.

Transparency

Declaration of financial/other relationships

ST, AS, KBW, KLT, and MHJ are current employees of and own shares in Eli Lilly and Company.

Author contributions

ST: study conception, study design, data acquisition, data analysis, data interpretation, drafting, and critical revision for important intellectual content. AS: data interpretation and critical revision for important intellectual content. KBW: study conception, data interpretation, and critical revision for important intellectual content. KLT: data acquisition, data interpretation, and critical revision for important intellectual content. MHJ: study conception, study design, data interpretation, drafting, and critical revision for important intellectual content.

Acknowledgements

The authors would like to thank Declan O’Dea and Sinead O'Loughlin from Eli Lilly and Company, and Keerthana Muthiah and Maanasa Surampally from Eli Lilly Services India Private Limited for their medical writing/editing support.

Reviewer disclosure statements

Peer reviewers on this manuscript have no relevant financial or other relationships to disclose.

Additional information

Funding

This work was sponsored by Eli Lilly and Company, Indianapolis, Indiana, USA.