Moving beyond the Cox proportional hazards model in survival data analysis: a cervical cancer study.

OBJECTIVES: This study explored the prognostic factors and developed a prediction model for Chinese-American (CA) cervical cancer (CC) patients. We compared two alternative models (the restricted mean survival time (RMST) model and the proportional baselines landmark supermodel (PBLS model, producing dynamic prediction)) versus the Cox proportional hazards model in the context of time-varying effects. SETTING AND DATA SOURCES: A total of 713 CA women with CC and available covariates (age at diagnosis, International Federation of Gynecology and Obstetrics (FIGO) stage, lymph node metastasis and radiation) from the Surveillance, Epidemiology and End Results database were included.
DESIGN: We applied the Cox proportional hazards model to analyse the all-cause mortality with the proportional hazards assumption. Additionally, we applied two alternative models to analyse covariates with time-varying effects. The performances of the models were compared using the C-index for discrimination and the shrinkage slope for calibration.
RESULTS: Older patients had a worse survival rate than younger patients. Advanced FIGO stage patients showed a relatively poor survival rate and low life expectancy. Lymph node metastasis was an unfavourable prognostic factor in our models. Age at diagnosis, FIGO stage and lymph node metastasis represented time-varying effects from the PBLS model. Additionally, radiation showed no impact on survival in any model. Dynamic prediction presented a better performance for 5-year dynamic death rates than did the Cox proportional hazards model.
CONCLUSIONS: With the time-varying effects, the RMST model was suggested to explore diagnosis factors, and the PBLS model was recommended to predict a patient's w-year dynamic death rate. © Author(s) (or their employer(s)) 2020. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

The present study explored the effect of variation in prognostic factors over time on overall survival and developed a prediction model to obtain the 5-year dynamic death rate for Chinese-American cervical cancer patients.

This is the first study to compare two alternative methods (the restricted mean survival time (RMST) model and the proportional baselines landmark supermodel (PBLS model, producing dynamic prediction)) versus the Cox proportional hazards model in the context of time-varying effects.

This study compared the performance of a variety of models using the C-index for discrimination and the shrinkage slope for calibration.

The use of the RMST model and PBLS model will provide statistical suggestions that could assess risks and predict dynamic death rate for clinicians and patients when the proportional hazards assumption fails.

Introduction

According to Surveillance, Epidemiology and End Results (SEER) Cancer Statistics Review, 13 170 estimated new cases of cervical cancer (CC, 0.7% of all new cancer cases) and 4250 estimated new deaths (0.7% of all new cancer deaths) will occur in 2019,1 making it the third most common cancer diagnosis and cause of death among gynaecological cancers in the USA. With the development of the human papillomavirus (HPV) vaccine, the incidence of CC will likely continue to decline both in the USA and worldwide.2 Unfortunately, despite these advances, Asian Americans are often reported to have low CC screening rates and high cancer mortality rates.3 4 Among the Asian-American subgroups, the Chinese-American (CA) subgroup is the largest population and the fifth fastest-growing racial population in the USA. However, few studies are available in the literature regarding exploration of the risk factors and prognosis of CC for CA women. Therefore, a survival model specifically for CA women with CC would be widely useful and clinically beneficial.

The Cox proportional hazards (PHs) model is the most common regression modelling framework to explore prognostic factors and to estimate survival rates. Importantly, the Cox PHs model has to satisfy the PHs assumption, which means that the ratio of the two hazard functions is constant over time. However, Trinquart et al 5 found evidence of the non-proportionality of hazards (ie, the HR is not constant over time) in 24% (13/54) of trials published in five leading oncology journals. In such situations, a constant HR (averaged HR over the follow-up) from the Cox PHs model may result in misleading conclusions.6 7

In recent years, except for the HR from the Cox PHs model, several authors5 8–11 have advocated the restricted mean survival time (RMST) difference as a robust and clinically interpretable summary measure of the survival time distribution that does not rely on the PHs assumption. By definition, the RMST is the mean survival time to a prespecified time point () and can be evaluated as the area under the survival curve S(t) up to this restricted survival time ().9 10 When the event of interest is death, it can be explained easily as ‘τ-year life expectation’. Thus, we also applied the RMST differences to explore prognostic factors in CA CC.

However, the above two models cannot be used to predict from an arbitrary prediction time point (s) during the follow-up period, as the models have been designed for using immediately after diagnosis.12 For example, patients may pay more attention to the survival probability or mortality of ‘w’ years after a cancer diagnosis, which is often asked by the question ‘How long will I live?’ or ‘What is the probability of being alive ‘w’ years from now?’. Furthermore, the question is asked at the start of the treatment as well as at any time during the follow-up visits. This question, which implies a conditional probability of surviving an additional w years, given that the patient is still alive at that time point,13 14 emphasises the dynamic use of the predictive model.15 Houwelingen and Putter16 17 designed the proportional baselines landmark supermodel (PBLS model; dynamic prediction) based on the Cox PHs model, considering time-varying effects. Thus, we can use this model to obtain individual death rates after an arbitrary time point s during follow-up. This concept of continually updating the w-year death rate from a certain s is referred to as the w-year dynamic death rate (DDR).

One aim of this article was to explore the effect of variation in prognostic factors over time on overall survival and to obtain the 5-year DDR for CA CC patients based on the SEER data. Another aim was to introduce two alternative methods when the PHs assumption failed and to provide statistical suggestions that could assess risks and predict survival probability for clinicians and patients.

Methods

Data sources

The SEER programme of the National Cancer Institute is an authoritative source of information on cancer incidence and survival in the USA that includes patient demographics, tumour characteristics, treatment utilisation and mortality. It is currently collecting and publishing cancer data from population-based cancer registries covering approximately 30% of the US population from 2013. The SEER database is an open public database, and the release of data from the SEER database does not require informed patient consent as patient privacy information is protected by the SEER cancer registries in every state of the USA. Access to the SEER 9 database was obtained after execution of A signed SEER Research Data Use Agreement (DUA).

Study design

Women diagnosed with CC were identified using the International Classification of Diseases for Oncology third Edition (ICD-O-3)18 codes, including cancer of the endocervix (C53.0), exocervix (C53.1), overlapping lesion of cervix uteri (C53.8) and unspecified cervix uteri cases (C53.9) for Chinese patients (863 cases) from 1988 to 2013. Patients with multiple primaries (0 case), autopsy only and death certificate only (4 cases), no survival time record (20 cases), not the first tumour (40 cases), unknown age at diagnosis (0 cases), unspecified International Federation of Gynecology and Obstetrics (FIGO) stage (36 cases), unknown lymph node metastasis (LNM) (75 cases) or unspecified radiation therapy (3 cases) were excluded (total number: 150 cases; some patients may miss several variables). Ultimately, 713 patients were involved in this analysis (see online supplementary eFigure 1).

Patient and public involvement

Patients were not involved in setting the research questions or planning the study. Researchers do not know the identities of the study participants.

Statistical analysis

The outcome of this analysis was all-cause mortality. The continuous variable (age at diagnosis) was indicated with median (range), and categorical data (FIGO stage, LNM and radiation therapy) were described using counts. The age at diagnosis was divided into per 10 years to better represent the covariate effect. The FIGO stage was dichotomised into categorical variables for easy statistical interpretation.

The Kaplan-Meier curves were graphed to describe overall survival. The univariable and multivariable Cox PHs models were used to estimate HRs and 95% CIs of the age at diagnosis, FIGO stage, LNM and radiation factors. Using the Grambsch-Therneau test,19 we found that the age at diagnosis and FIGO stage violated the PHs assumption, in which cases the interpretation of a constant HR was difficult. Therefore, RMST was used as another appropriate outcome measure in our analysis. We chose a restricted survival time () first and then computed the pseudo-value for each observation and finally fitted a generalised linear model (GLM) to assess the prognostic factors.20

For dynamic prediction analysis,16 we first fixed the prediction window at w=5 years—that is, the ‘w’ in response to the patient’s question ‘How long will I live ?’ at any prediction time point (s). Next, we selected a set of prediction landmark time points {s 1…,s 101} from prediction time point s at approximately every 2.4 months between 0 and 20 years after the diagnosis of CC. For each landmark time point s, we constructed a landmark dataset R, which consisted of individuals still alive at s, with administrative censoring at s (see online supplementary eFigure 2). These landmark data sets R were stacked to create a super prediction data set R. We then applied a PBLS model to the stacked super dataset R, which allows the baseline hazard to change smoothly with the landmark time point. This model can be fitted directly by applying a simple Cox model to the dataset R, in which the landmark time point is used as a covariate. Using this analysis, we can obtain the 5-year DDR at any prediction time point, s, between 0 and 20 years after diagnosis of CC. In addition, to test for time-varying covariate effects, interactions between covariates and s (both linear and quadratic) were also included in the PBLS model, and the Akaike information criterion (AIC) was used to decide the final model for dynamic prediction.

The performances of the prediction models (the Cox PHs model and the PBLS model (w=5)) were evaluated in terms of both discrimination and calibration. The model’s capacity to correctly discriminate among patients was evaluated using Harrell’s C-index, while the calibration using the heuristic shrinkage factor was based on 10 000 bootstrap resamplings.

All analyses employed R 3.2.4 software, and the significance level was set to 0.05.

Results

Characteristics of the study population

The characteristics of the 713 CC patients with 162 (22.7%) deaths, with a long-term (nearly 26 years, ie, the maximum follow-up time in this data) survival rate of 69.4% (95% CI 64.5% to 74.7%), are summarised in table 1. The median of age at diagnosis was 42 (15–97) years. The median death time and median follow-up time were 4.08 (range: 0.08–24.08) and 16.67 (range: 0.08–25.92) years, respectively.

Table 1

Two models for overall survival

Variable	Baseline characteristics	Cox proportional hazards model			RMST based on pseudo-values (GLM)
	No. of patients (deaths)	Univariable HR (95% CI)	Multivariable HR (95% CI) *	P value†	Univariable RMST(95% CI)‡	Multivariable RMST(95% CI)*	P value†
		–	–		‡	22.12 (20.18 to 24.05)	<0.001
Age at diagnosis (per 10 years)	713 (162)	1.90 (1.72 to 2.10)§	1.59 (1.42 to 1.79)¶	<0.001	−1.79 (−2.11 to –1.47)	−1.06 (−1.43 to –0.69)	<0.001
FIGO stage
In situ	375 (44)	0.29 (0.19 to 0.44)§	0.56 (0.35 to 0.88)¶	0.013	3.01 (2.01 to 4.02)	1.35 (0.30 to 2.41)	0.012
I (reference)	205 (46)
II	101 (50)	3.28 (2.19 to 4.92)	1.57 (0.97 to 2.54)	0.067	−5.07 (−6.98 to –3.17)	−3.04 (−5.34 to –0.73)	0.010
III	32 (22)	7.23 (4.32 to 12.12)	2.26 (1.23 to 4.15)¶	0.009	−9.64 (-12.70 to –6.58)	−6.21 (−9.50 to –2.92)	<0.001
Lymph node metastasis
No (reference)	659 (132)
Yes	54 (30)	6.07 (4.04 to 9.12)	2.52 (1.61 to 3.94)	<0.001	−7.93 (-10.31 to –5.55)	−3.63 (−6.29 to –0.97)	0.008
Radiation
No (reference)	532 (77)
Yes	181 (85)	7.36 (5.32 to 10.19)	1.43 (0.87 to 2.33)	0.156	−6.93 (−8.25 to –5.61)	−1.09 (−3.05 to 0.86)	0.273

*Adjusted for the age at diagnosis, FIGO stage, lymph node metastasis and radiation.

†P values were obtained from the Wald test and adjusted for the age at diagnosis, FIGO stage, lymph node metastasis and radiation.

‡The intercepts for each variable are 24.96, 16.23, 17.26 and 18.42, respectively.

§Grambsch-Therneau proportional hazards test, P<0.05 for univariate analysis.

¶Grambsch-Therneau proportional hazards test, P<0.05 adjusted for age at diagnosis, FIGO stage, lymph node metastasis and radiation.

FIGO, International Federation of Gynecology and Obstetrics; GLM, generalised linear model; RMST, restricted mean survival time.

Effects of prognostic factors on overall survival

According to the Cox PHs model, all covariates showed significant differences (P<0.001) in the univariate analysis, while the multivariate analysis demonstrated that radiation did not have a significant effect (P=0.156) on CC patients (HR=1.43; 95% CI 0.87 to 2.33) (for more details, see online supplementary eFigure 3 and table 1). However, the Grambsch-Therneau test showed that the age at diagnosis and FIGO stage were covariates with time-varying effects, indicating that constant HRs for these covariates obtained from the Cox PHs model were unreliable.

Table 1 also shows the RMST differences. The patients without LNM had a longer restricted mean survival time, namely, 3.63 years, than the women presenting with LNM after 20 years of follow-up in the multivariate analysis. The women with FIGO stage II compared with patients with stage I showed a statistically significant RMST difference (−3.04; 95% CI −5.34 to −0.73), while the HR (1.57; 95% CI 0.97 to 2.54) showed no significance in the multivariate analysis.

Table 2 demonstrates the regression coefficients in the PBLS model, and figure 1 shows the HR 5(s)s, that is, the HR from any prediction time point s to s+5. The age at diagnosis revealed a changing impact on HR 5(s) with each successive prediction time point (s), indicating that CC women had increased 5-year death rates over time (figure 1A). The women in FIGO stage I compared with stage in situ illustrated a significant time-varying effect, while the stage II versus stage I and stage III versus stage I did not (figure 1B1–B3); however, the overall time-varying effect was significant. The LNM also had the potential to change over time in the 5-year dynamic HR, HR 5(s), which could be calculated by the following formula:

Table 2

The dynamic prediction used the PBLS model (w=5 years)

	Variable	Time function*	Coefficient	SE	P value†	P value‡
β LM (s)	Age at diagnosis(per 10 years)	1	0.18	0.11	0.082	<0.001
		s/20	1.14	0.25	<0.001
	FIGO stage(1994 edition)					<0.001
	In situ	1	−2.39	0.55	<0.001
		s/20	11.10	2.50	<0.001
		(s/20)2	−9.07	2.56	<0.001
	I (reference)
	II	1	0.39	0.40	0.335
		s/20	2.68	2.79	0.336
		(s/20)2	−3.44	3.14	0.274
	III	1	1.08	0.44	0.015
		s/20	2.20	4.69	0.639
		(s/20)2	−18.40	11.22	0.101
	Lymph node metastasis					<0.001
	No (reference)
	Yes	1	0.40	0.40	0.316
		s/20	5.66	3.15	0.072
		(s/20)2	−9.56	4.00	0.017
	Radiation
	No (reference)
	Yes	1	0.34	0.30	0.260	0.252
θ(s)	s/20	−10.60	1.89	<0.001
	(s/20)2	3.50	1.55	0.024

* Time function: ,

† P value for the time-varying trend.

‡ P value for the whole time-varying effect adjusted for age at diagnosis, FIGO stage, lymph node metastasis and radiation.

FIGO, International Federation of Gynecology and Obstetrics; PBLS, proportional baselines landmark super.

Figure 1

Dynamic HRs with 95% CIs in the dynamic prediction PBLS model (w=5 years). The red line represents HR=1. The black solid line represents dynamic HRs, and the black dashed line represents 95% CIs. FIGO, International Federation of Gynecology and Obstetrics; PBLS, proportional baselines landmark super.

The HR 5 (s) (figure 1C) increased from the start of the diagnosis but decreased 6 years after the diagnosis and was nearly 0 with a 20-year history of CC (, , , ). Radiation was not statistically significant (figure 1D).

Individual dynamic prediction

The 5-year DDRs for patients with different baseline characteristics (see online supplementary eTable 1) are shown in figure 2. For example, the green line in figure 2 represents the probability of dying within 5 years for patient B with the following characteristics: diagnosed at the age of 45 years, FIGO stage I, without LNM and no radiation. If still alive at 7 years after diagnosis, that is, the landmark time equaled 7 years, the probability of dying within 5 years was less than 0.05 for patient B.

Figure 2

Changes in the 5-year dynamic death rate in nine example patients. FIGO, International Federation of Gynecology and Obstetrics.

Model assessment

The Cox PHs model and the PBLS model (w=5) presented different discrimination capacities (0.836 and 0.868, respectively, for the C-index). The PBLS model demonstrated a better discriminative capacity than the Cox PHs model. The calibration slopes (0.958 and 0.996, respectively) for these two models were similar for the apparent performance, while the PBLS model showed a better fit.

Discussion

We used the Cox PHs model, as well as two alternative methods (the RMST difference model and the PBLS model), to explore the prognostic factors and provide 5-year DDR predictions for CA CC patients. To our knowledge, these two alternative models for CA CC patients have not been previously reported.

Regarding factor exploration, the Cox PHs model revealed that age at diagnosis, FIGO stage and LNM were independent risk factors of the prognosis of CA CC patients. However, the presence of non-PHs indicated that the effects of some covariates were changeable over time. In such a situation, RMST, the area under the survival curve S(t) up to the restricted time point (τ), has great potential as another meaningful and sensitive outcome measure in the analysis with a time-to-event outcome. Compared with HR, the RMST difference does not rely on the PHs assumption and is readily interpretable. For instance, women without LNM had a longer mean survival time than the patients with LNM after the 20-year follow-up (blue area in online supplementary eFigure 4), namely, 7.93 years in the univariate analysis. In addition, the RMST is estimable even under heavy censoring, while the CI for HR from the Cox PHs model can be unstable and can vary widely. However, satisfactory results of RMST analysis must depend on the choice of the restricted time point (τ). Generally, the value of τ should be prespecified to avoid selection bias (after seeing the data). An appropriate time point should consider both the statistical method and the clinical reality. For example, in a trial on a severe cancer, 1-year or 2-year survival may be an appropriate measure, so τ=1 year or 2 years would be reasonable, while for cancers with a relatively longer clinical course, such as CC, τ=10 or even 20 years is more realistic. In the absence of such a choice, a default τ may be taken to be slightly below the maximum expected follow-up time.9 In summary, the HR and RMST difference are complementary techniques that provide alternative methods of summarising treatment effects.8

There are also some other models that are commonly used when the PHs assumption fails, such as the accelerated failure time (AFT) model and the transformation model. However, these two models have certain limitations in their use and their results are not easy to interpret. For example, although the AFT model is simple and does not need to meet the PHs assumption, the distribution of survival time must be specified in advance, and random errors, noise and covariates are assumed to be independent and so on. Due to the special distribution of the survival time in the data in this study (for more details, see online supplementary eFile 1), the results may not be accurate. The transformation model also has some advantages and has proved its relevance in the two special cases (the Cox PHs model and the proportional odds model; for details, see online supplementary eFile 1), but for other choices of r, the regression coefficients are more difficult to interpret because they refer to the scale given by the unknown h. Therefore, the two models are introduced in online supplementary eFile 1, and their results are shown for reference only (see online supplementary eTables 2 and 3, eFigure 5).

Exploring the prognostic factors is important for clinical research, but a good prognostic prediction model might be more important to patients. We used the PBLS model to perform dynamic prediction and it had several strengths over the above two models. First, the PBLS model can effectively find the time-varying effect of covariates for both estimation and prediction (see table 2 and figure 1). Second, patients may focus more on ‘What is the probability of still being alive w years from now?’ during the follow-up. The PBLS model can provide the HR (s) or w-year DDR at different prediction times. The Cox PHs model in this analysis could only provide a static HR for the whole follow-up and obtain the 5-year survival rate from the diagnosis time to 5 years. Furthermore, the discrimination and calibration in the PBLS model were clearly larger than those in the Cox PHs model, indicating that PBLS is more suitable in reality than the Cox PHs model for prediction when the PHs assumption is unsatisfied. Third, dynamic prediction is based on landmark analysis; thus, this model could overcome the guarantee time bias,21 especially in the group classification.

However, dynamic prediction must be considered cautiously. We noticed that the dramatic 5-year dynamic HR for the dummy variable (FIGO III vs FIGO I) was less than 1 after 7.4 years. The reason is that, at each landmark time point (s), we have a corresponding dataset (R) (see online supplementary eFigure 6). As the landmark time point increased, the patients whose survival times were less than s would be excluded, leading to a decrease in the number of patients in the dataset R. The cumulative distribution function (F(s)) (see online supplementary eFigure 7) was calculated by dividing the number of deaths by the number of patients in each dataset. It can be seen that, for FIGO III, F(s) became 0 after a 7.4-year follow-up as the number of deaths equals 0. Thus, the reason for the dramatic HR 5(s) less than 1 after 7.4 years was apparent. Therefore, we could not simply state that the patients in FIGO III had better survival than those in FIGO I at the late follow-up. The same situation occurred regarding LNM after the 11-year follow-up.

In addition, several limitations of the current study must be noted. First, this is a retrospective observational study, and consequently, some bias may have affected the analysis. Second, SEER database lacks some important variables related to patient status, such as HPV infection, the information of chemotherapy, and economic conditions, which are associated with prognosis. Third, because SEER does not include information on cancer recurrence, we are unable to analyse other endpoints, such as time to recurrence, or disease-free survival. However, because cause-of-death information is available, an additional cause-specific analysis was performed to model CC mortality (for details, see online supplementary eFile 2), and the results are shown in online supplementary eTables 4–7, eFigures 8 and 9.

Above all, older patients had a worse survival rate than younger ones, a finding similar to that of Nghiem VT for Chinese patients. The FIGO stage is an essential prognostic factor in CC patients, and our study confirmed that advanced FIGO stage patients had a poor survival rate.22 23 LNM was an unfavourable prognostic factor in our models. Several studies24 25 have indicated that patients with LNM have a rather poor prognosis. The new finding that age at diagnosis, FIGO stage and LNM have time-varying effects was found from the dynamic prediction. In this study, radiation had no impact on survival in CA CC, a finding that differs by race in some reports.26 27

In summary, the Cox PHs model is unsuitable for survival analysis when the PHs assumption is unsatisfied; thus, the RMST difference as another effective outcome measure that overcomes the PHs handicap, which would help clinicians to explore the effects of prognostic factors for CA CC, and the PBLS model (dynamic prediction) permits continuous revision of a CA CC patient’s w-year DDR at the start of the diagnosis as well as at different prediction times s during follow-up. We also recommend the latter two alternatives considering time-varying factors when the PHs assumption is not satisfied. But in actual use process, we must pay attention to the choice of the restricted time point (τ) in the RMST model and the interpretation of the PBLS model.