Predictive Modeling of Long-Term Survivors with Stage IV Breast Cancer Using the SEER-Medicare Dataset

1. Introduction

A total of 43,170 US women died from breast cancer (BC) in 2023, most from stage IV disease^1-5. Despite some progress in prolonging survival, their overall median survival remains 2-3 years^3,4. Treating patients with stage IV BC can prolong population averaged survival by only a few months^6-8. Systemic treatment options vary in the wide range of scenarios generated by patient- and cancer-related variables, and include many possible options from combinations of chemo-, hormone-, bio-, or immune-therapy categories^{4, 9,10}.

However, despite the grim overall numbers, there are some long-term survivors^11,12. There are several variables linked to survival times. Key among these is the presentation of the disease. Patients with de novo presentation of stage IV BC have a median survival superior to that of patient with recurrent disease^13-19. In one study, the median survival of patients with de novo metastatic BC was 41 months while the median survival of patients with recurrent BC was 25 months¹⁷. In patients with recurrent disease, shorter metastasis-free intervals of <24 months are associated with worse survival than longer metastasis-free intervals^14,16,18,19. Neoadjuvant treatment for the original local disease had a significant positive association with survival²⁰.

In patients with de novo diagnosis of stage IV disease, risk stratification models identified variables at the time of diagnosis positively and negativley associated with overall survival. These variables included age at diagnosis, particularly age over 80 years, premenopausal status, Karnofsky performance status, estrogen receptor (ER)-, progesterone receptor (PR)- and human epidermal growth factor receptor 2 (Her2)-status, tumor histology, grade, tumor size categories 0-3 vs. 4, lymph node status, number of metastatic organ systems from 1-4, number of metastases, number of metastatic sites, bone only metastases +/-, visceral metastases +/-, brain only metastases +/-, frontline therapy with Her2 antibody or hormone blockade, and combinations of these variables^14,17,21-28.

In addition to these time-fixed variables, time-varying variables occurring after initial diagnosis of de-novo stage IV disease are predictive of survival. These include resection of the primary tumor,^29-32, of locoregional metastases³³, or of lung only metastases³⁴. The impact of primary tumor resection on survival is affected by age, race, cohabitation status, income, tumor size, grade and histology, ER, PR or Her2 status, expression of Ki₆₇ and CA₁₅, alkaline phosphatase, lymphovascular invasion, lymph node status, metastasis to brain, liver, and lung, and chemotherapy^27.29,30.

With respect to systemic therapy, although population- averaged survival extension is modest, there are multiple treatment- and circumstance-associated combinations that can result in prolonged survival in specific cases. Systemic treatment options consist of individual or combinations of categories of chemo-, hormone-, bio-, or immune-therapy^4,9,10. Variables that impact survival include treatment timing, dosing intensity and density, duration, sequential or intermittent drug administration, alternating modalities, first or subsequent lines of treatment^11,12, adverse events (AEs)⁹, particularly in association with other negative patient risk factors^35-42. For exampler, long-term survival of Her2+ patients with Her2-targeted therapy, taxane-based, hormone maintenance therapy, nab-paclitaxel therapy, or multimodality therapy, with or without local treatment was dependent on hormone receptor expression, disease burden, soft tissue or bone vs visceral metastases, resection of the primary tumor or metastases and young age^43-49, and on freedom from disease progression after a year of therapy⁵⁰.

Overall, the highest responses in stage IV BC are observed with the first line of therapy, followed by rapidly diminishing returns^11,12,42. Disease progression after the first few treatment regimens results in exceptionally heterogeneous disease, with cancer cells and tumor stroma selected for treatment resistance and an enhanced capacity to re-metastasize⁵¹. Decisions to treat to prolong survival also include considerations of treatment-induced AEs^52-61. Therefore, it would be highly beneficial to patients and their caregivers to know the patient- and cancer-spcecific circumstances where a patient could expect long term survival benefits from specific treatmets. Population-based predictive normograms achieved concordance indices (C-index) below 0.700 using significant univariate risk factors^17,62, which rose to 0.737 using Cox-regression modeling of five prognostic categories⁶³. However, predictive models for individual patient survival are needed to navigate the exceptionally large number of potential circumstance specific scenarios.

Given the rich and variable population-based predictive data with significant variation and the general limited survival with stage IV BC, we have undertaken a Deep Learning (a subfield of artificial intelligence) approach to generate predictive modeling of individual patient survival from the SEER-Medicare data of time-fixed and time-varying covariates. Using this unique approach, we can identify patient and treatment scenarios that will result in long term survival. The models will help caregivers identify specific treatments from a wide variety of choices that will optimize survival in individual patients and unique circumstances.

2. Methods

2.1. Study Data: SEER-Medicare Linked Dataset

Unlike the SEER dataset, which does not provide data on cancer therapy, the SEER-Medicare (S-M) linked dataset created almost 30 years ago is considered a major source to assess cancer care and outcomes in the US⁶⁴. The S-M linked dataset provides information about Medicare beneficiaries with cancer. This dataset includes cancer incidence data for about 26% of the US population in various regions. Medicare data have SEER cancer file, a person-level file that provides SEER demographic and clinical information for up to 10 primary cancer diagnoses, treatments, and mortality. Medicare files capture the fee-for-service claims from hospitals (MEDPAR), outpatient facilities (Outpatient), National Claims History (NCH), hospice care, home health agencies, and the Master Beneficiary Summary File (MBSF), which contains Medicare enrollment information. The CCflag file includes, for every patient, the date the patient was diagnosed with one of 22 chronic conditions. We use these data to compute a time-varying patient's comorbidity index. There were 883,053 BC patients during 1991-2016. The number of cancer patients and their corresponding records in each file were outlined previously⁶⁵. The patient numbers in each file reflect the fact that some patients had recorded encounters in multiple files. Based on the Observational Medical Outcomes Partnership (OMOP) common data model, we fused the data in the various files. This enabled us to transform data contained within those disparate files into a common format (data model) as well as a common representation (terminologies, vocabularies, coding schemes).

2.2. Inclusion Criteria

In our study, we included women with the diagnosis of stage IV breast cancer, who have not had any other malignancy history except non-melanoma skin and eyelid cancer, as a standard in NCI clinical trials⁶⁶, included all comorbidities recorded at every visit from the date of diagnosis, averaged over the course of all the treatments assessed for each patient, all prior treatments, delineated age, race, marital status, breast cancer stage, grade, hormonal and Her2 status, and laterality. We included only patients who were enrolled in both Parts A and B of Medicare with no HMO enrollment 1 month prior to diagnosis with breast cancer and whose age at enrollment was 65 or older who qualified by age and not disability. The enrolled population was, therefore, skewed by age and only represented an elderly subset of the breast cancer patient population. We included patients who were enrolled in both Parts A and B of Medicare with no HMO enrollment from 1 month prior to diagnosis through 20 years following diagnosis, hospice or death to ensure that subjects were continuously enrolled in the proper parts of Medicare during the study period, as before⁶⁷. The number of entries, patients, mean age and comorbidity index are outlined in Table 1.

2.3. Data Cleaning and Standardization

Data cleaning was implemented to remove duplicate records in the dataset. We included patients who have valid patient id, cancer type (breast, or breast plus skin, or eyelid), stage (I-III, IV), sex, date of diagnosis (month/year; in case of a missing day of the month, we assumed the first day of the month). A patient visit (inpatient, outpatient, and carrier claims) was included if it had a valid patient id and a valid date. We included only valid diagnoses, procedures, or HCPCS codes at a given visit. We applied the Melt transformation to convert a visit record (containing a mix of valid and missing codes) from wide format to stacked/long form, thus being able to delete invalid codes while keeping only valid ones for each visit. We performed data transformation and standardization, ensuring all features were numeric and standardized, and applied bucketing for continuous values, e.g., age and comorbidity index.

2.4. Treatments and Adverse Events: Encoding and Embeddings

We included all patients' ICD-9 and HCPCS codes. We annotated 141 HCPCS drug J codes: 82 chemotherapy drugs, 49 biotherapy drugs, and ten hormone therapy drugs, and consolidated related adverse events into 18 categories to enable a more comprehensive analysis of the variables associated with their occurrence, as before ⁶⁷.

2.5. Predictive Modeling

2.5.1. Discrete Time-to-Event Data

In our dataset, as is in real life, patient follow-up visits occur on a given day with irregular gaps between two consecutive visits. A patient has two sets of covariates: time-fixed covariates (e.g., age at diagnosis) and time-varying covariates (e.g., current age, current comorbidity index, treatments administered at this visit and earlier visits, and adverse events). A patient's survival status is recorded at each visit while the patient is at risk, i.e., has not yet experienced the event (death in our case) of interest. We have 14,312 breast cancer stage IV patients and 1,880,153 entries over the period 1991-2015. Since we selected patients with ICD-9 diagnosis codes for this study, the study period was set to 1991-2015. Thus, our study is right-censored, i.e., some patients are not followed all the way to their event (death) time, resulting in censored time (2015-12-31) instead of event time. Thus, instead of observing the event (death in our case) at time, $T^{*}$, we observe a possibly right-censored event time, $T=\min \left\{T^{*}, C^{*}\right\}, where C^{*}$is the censoring time (2015-12-31 in our case). We also observe the indicator $D=1\{T=T^{*}\}$ labeling the observed event time as an event or a censored observation.

The general idea of prediction using this discrete-time framework is to build models predicting the likelihood of surviving each visit. In our case, the event of interest is death. Let $T^{*}$ be the time to death since diagnosis, and $T=\min \left\{T^{*}, C^{*}\right\}, where C^{*}$is the censoring time, and $\mathit{x}$ is a covariate vector. We are predicting the probability of a patient experiencing failure (death) by time $t$. This probability, referred to as the cumulative incidence function, is given by $P\left(T\le t|\mathit{x}\right).$ An alternative to this probability is the survival function, $S\left(t\right)=P\left(T>t|\mathit{x}\right)$ and the hazard rate, $h\left(t|\mathit{x}\right).$The Cox proportional hazards (CPH)⁶⁸ is the most used Model in survival analysis. It provides a semi-parametric hazard rate, $h\left(t|\mathit{x}\right)$ that is the product of a baseline hazard, $h_0\left(t\right)$ and a relative risk function, $e^{g\left(\mathit{x}\right)}$, i.e., given by: $h\left(t|\mathit{x}\right)=h_0\left(t\right)e^{g\left(\mathit{x}\right)}, g\left(\mathit{x}\right)={\beta }^T\mathit{x}$, where $\mathit{x}$ is a covariate vector, and $\beta$ is a parameter vector. The Model's proportionality hazard assumption, i.e., the effect of each patient covariate is the same at all values of the follow-up time, is unrealistic for most clinical situations⁶⁹.

We applied the existing DL-based models DeepSurv⁷⁰ to extend the Cox regression models with neural networks estimators^71,72, and DeepHit^71-74, Nnet-survival⁶⁹, and Cox-Time⁷⁴ to extend the Cox regression model with neural networks and overcome the proportionality assumption of the Cox model, as before⁷⁵.

2.5.2. Time-Varying Covariates

The above methods, while dealing with discerete-time data they assume the covariates$, {\mathit{x}}_{\mathit{i}}$ of a patient, $i$, are time-fixed. In our case, a realistic prediction of the patient's survival needs to take into account not only the patient's marital status, race, hispanic, age, laterality, grade, hormonal and Her2 status, and laterality at the time of diagnosis, but also, the administered treatments, induced adverse events, and comorbidity index, and the age at each visit. Thus, being able to predict associations between specific treatments in specific scenarios and outcomes relating to quality of life, progression-free survival, and survival in the spectrum of stage IV cancer settings.

To achieve this objective, we propose to extend a patient covariate vector$, {\mathit{x}}_{\mathit{i}}$ to include not only time-fixed covariates but also covariates that summarize the patient's history from previous visits. Specifically, for a given treatment$, {TR}_{\mathit{j}}for j=1, 2,\dots , 46$, ${TR}_{\mathit{i}\mathit{j}}$is a tally of the number of times ${TR}_{\mathit{j}}$ was administered to this patient, $i$, from the time of the diagnosis to the time of death/end of the study. Similarly, for a given adverse event inuced$, {AE}_{\mathit{k}}for k=1, 2,\dots , 18$, ${AE}_{\mathit{i}\mathit{k}}$is a tally of the number of times patient, $i$, experienced ${AE}_{\mathit{k}}$, from the time of the diagnosis to the time of death/end of the study. We divide the age of a patient into 6 bins: $bi{n}_1$<=65, 65< $bi{n}_2$<=70, 70< $bi{n}_3$<=75, 75< $bi{n}_4$ <=80, 80< $bi{n}_5$ <=85, $bi{n}_6$ >85. We have then, ${AGE}_{\mathit{i}\mathit{b}}$is a tally of the number of times the age of patient, $i$, falls within the $bi{n}_b for b=1, 2, \dots , 6$ , from the time of the diagnosis to the time of death/end of the study. We handle a patient comorbidity index similar to a patient age. We divide the comorbidity index of a patient into 6 bins: $bi{n}_1$<=2, 2< $bi{n}_2$<=4, 4< $bi{n}_3$<=6, 6< $bi{n}_4$ <=8, 8< $bi{n}_5$ <=10, $bi{n}_6$ >10. We have then, $COM{B}_{ib}$is a tally of the number of times the comorbidity index of patient, $i$, falls within the 6 comorbidity bins: $1, 2, \dots , 6$ from the time of the diagnosis to the time of death/end of the study.

2.6. Experiments

2.6.1. Data Prepration

Our data consisted of patient-level records with a single row for each patient and a column for each patient covariate as described above. A patient had a column indicating the corresponding value of $T=\min \left\{T^{*}, C^{*}\right\},$and a column indicating the patient event status at that time. The time horizon (1991-01-01 through 2015-012-31) was divided into months.

Our Model provides an estimate of a patient's complete survival curve, S(t); these estimated survival curves represent event probabilities as time functions.

2.6.2. Performance Metrics

In general, metrics for measuring the Model's predictive performance are often in terms of discrimination, i.e., the Model's ability to distinguish between patients with high and low risk of experiencing the event, and calibration, i.e., the agreement between the estimated and actual incidence of the event⁷⁶.

At time 0, patients are event-free, and their survival status changes from one time period to another. Thus, performance metrics for assessing the predictive accuracy of time-to-event outcomes are time-dependent. In addition to accounting for time, performance metrics must also account for censoring. Two such metrics that are commonly used to measure the predictive accuracy of survival prediction models are: 1) time-dependent concordance index⁷⁷, and time-dependent Brier score Brier score^78-80. The concordance index, or C-index, referred to as the AUC (area under the receiver operating characteristics curve), is the most commonly used method to measure the Model's predictive accuracy. It estimates the probability that, for a randomly selected two patients, the predicted survival times of the two patients have the same ordering as their true survival times⁸¹. The C-index is a scale-free measure with 1 representing perfect discrimination and 0.5 representing discriminative ability similar to chance. Brier score, the predictive error, assesses the Model's calibration and discrimination ability, with lower values indicating better prediction. Similar to⁷⁴, we use the time-dependent C-index by Antolini et al.⁷⁷ and will use the integrated Brier score (IBS) by Graf et al.⁷⁹ distribution.

2.6.3. Hyperparameter Tuning

Hyperparameters are neural network parameters fixed by design and which accounts for censored patients by weighting the score by the inverse of the estimated censoring not tuned by training but should be optimized. Hyperparameter optimization is an important part of deep learning, in general. We applied Amazon SageMaker Python SDK (software development kit), an open-source library, to fine-tune the model by identifying optimal values of the network's hyperparameters, as before⁷⁵. Table 2 includes the list of the hyperparameters and their ranges. The experiments will be conducted using five-fold cross-validation using a Bayesian optimization search scheme⁸². The Bayesian optimization scheme has been shown to outperform other state of the art global optimization algorithms on a number of challenging optimization benchmark functions⁸³. Figure 1 depicts the performance metric (time-dependent concordance) and each of the model's hyperparameters: n-layers, Ir, batch size, epochs, dropout, n_nodes, alpha, sigma.

• Validation

The base model (vs. the extended Model)

We used the pycox package⁷⁴. To validate our implementation, we applied the four models: DeepSurv, DeepHit, Nnet-survival, and Cox-Time to the real datasets METABRIC⁸⁴ and SUPPORT⁸⁵, which is publicly available on the Vanderbilt Biostatistics Web site. The experiments were conducted by five-fold cross-validation. As we demonstrated before⁷⁵, the published results of each of these models using the Metabric and SUPPORT datasets are within our results' 95% confidence level, respectively (we used only the time-dependent C-index as the performance metric)⁸⁶.

3. Results

3.1. Improved Predictive Capacity When Adding Time-Variable Covariates to Time Fixed Covariates

The results are presented in Table 3. In terms of concordance, we observe the significant improvement of the model with the proposed extended patients' covariates compared with that of the patients' time-fixed covariates. For example, the DeepSurv model's prediction error is 4% when using the proposed extended patients' covariates versus over 32% when using the patients' time-fixed covariates.

3.2. Predicited Survival Varies by Patient Covariates

The median survival of the five hypothetical patients ranges from less than one month to 100 months. The predictive models utilized the unique circumstances created by the interaction of the individual patient, cancer and treatment variables to generate the individual patient’s predicted survival curves. We cannot provide patient details because our DUA with the NCI SEER-Medicare does not permit including specific patient identifying information.

4. Discussion

Our results represent two significant achievements in predictive modeling of patient outcomes. The first achievement presents a unique and compelling opportunity to improve the prediction performance of the four deep learning models DeepSurv, DeepHit, Nnet-survival, and CoxTime, which outperform the Cox Proportional Hazard Model in survival analysis. Our models handle discrete-time distribution by extending the patients' covariates vectors to include both time-fixed covariates and covariates that summarize the patient's history from previous visits. This lowers the prediction error rate from 28-38% to 2-12% using the four deep learning models. The IBS can be viewed as the mean square error of prediction; lower values of the IBS indicate better predictive performance. The cause-specific time-dependent Concodrance-index accounts explicitly for censoring, and estimates the model's prediction error.

The second achievement is the application of the models to generate individual patients’ predictive survival curves, based on their unique patient, cancer and treatment features. This will permit predictive modeling in the myriad scenarios encountered with individual patients on likelihood of survival, AEs and impact of progressive therapy. The event pattern-based representation of associations relationships will enable us to perform reasoning and what-if analysis using temporal logic-based approaches; e.g., treating a patient who has progressed after multiple lines of systemic therapy and selecting the optimum treatment type most likely to prolong survival in the specific scenario generated from the conglomeration of the unique patient, cancer and treatment history variables. Outcomes from these models include the first demonstration that the occurrence of certain AE categories have a negative impact on survival^65[Ada24]. In that light, the models will inform patients and treating physicians of the most efficacious next line of treatment to prolong survival, or of the likelihood of experiencing primarily adverse events with a minimal likelihood of prolonging survival. The concordance index and the Integrated Brier scores validates the accuracy of individual patient predicted survival curves. Prospective studies integrating clinical and genomic data may identify unique clinicogenomic features of MBC patients who can achieve durable disease control without prolonged chemotherapy^47[Shi23].

There are several potential limitations to this study, however. Our patient population consisted of Medicare-enrolled patients and is restricted to patients over 65 years who qualified by age and not other medical conditions. Their age range, therefore, is not representative of the general population and their overall expected survival may be less than patients with stage IV BC of a younger age. We have consolidated AEs and treatment codes into general related categories to enable analysis. Thus, the impact of treatments and adverse events we have analyzed may vary among individual drugs or AEs collated into single categories. The AEs we considered as time-varying events are not graded in the SEER-Medicare data as minor or moderate vs. severe. Our future investigations based on data-driven hypotheses generated from this work will include a classification of non-severe and severe AEs, with the latter assigned to AEs associated with emergency department or hospital admissions or death.

5. Conclusions

These investigations demonstrate a significant improvement of the deep learning models DeepSurv, DeepHit, Nnet-survival, and CoxTime. These models outperform the Cox Proportional Hazard Model in survival analysis with the proposed extended patients' covariates compared with that of the patients' time-fixed covariates, reducing the error rate of concordance indices from 28-38% to 2-12%. The project also developed individual patients’ predictive survival curves, based on their unique patient, cancer and treatment and adverse events variables. The models provide a unique tool to medical care givers to generate realistic survival probabilities for individual patients with stage IV BC given their unique, specific circumstances. It will be a highly useful adjunct to the clinical decision-making in individual patients guiding potential next line of therapy.