Investigating Time-Varying Predictor Effects on Cardiovascular Outcomes in Breast Cancer Survivors.

INTRODUCTION
Breast cancer (BC) survival rates have improved over the years because of earlier detection and advances in treatment [1]. Following improved BC survival, several studies have reported an increase in cardiovascular disease (CVD) events, with a higher incidence observed among patients with BC than among those without BC [234]. The underlying reasons for this association are multifactorial. Certain BC treatments, such as cyclophosphamide, epirubicin, and trastuzumab, are known to be cardiotoxic [56], and all are commonly used to treat breast cancer. The European Society of Cardiology has suggested that identifying patients at a high-risk of CVD before initiating cancer therapy could allow clinicians to tailor treatment strategies and thereby reduce the risk of cardiotoxicity [7]. Although this has been examined in other studies, the extent to which the relevant factors exhibit a shifting prognostic impact over time has received far less attention. By applying Cox and alternative survival model specifications, in addition to a random survival forest (RSF) analysis, we demonstrated how pronounced time-varying effect can be when discussing time to CVD events among BC patients. Studies such as Boekel et al. [8] have described the overall incidence of CVD among patients with BC. However, the incidence and underlying risk factors may vary significantly across specific CVD types, such as atrial fibrillation compared with angina pectoris [9], underscoring the need for a more nuanced investigation. Therefore, we investigated how subtype-specific CVD outcomes are associated with prognostic factors in patients diagnosed with BC.

METHODS

Study endpoints
We examined the time to the first occurrence of six CVD endpoints: myocardial infarction (MCI), atrial fibrillation (AF), embolic stroke (ES), angina pectoris (AP), embolism, peripheral vascular disease (EPV), and heart failure (HF), with a particular focus on their associations with molecular subtype. The time to the first occurrence of MCI is hereafter referred to as the time to MCI, similar to the other CVD endpoints. Among the six CVD endpoints, HF represented the final and most clinically significant CVD in terms of cardiotoxicity. In addition, time to HF was considered the most relevant endpoint owing to its time-varying association with the patients’ molecular subtype. Hence, our primary priority was the HF endpoint.

Molecular subtype groups and BC treatments
Patients with BC were classified into four molecular subtype groups: 1) hormone receptor (HR)-positive and human epidermal growth factor receptor 2 (HER2)-negative (HR+/HER2−), 2) HR-positive and HER2-positive (HR+/HER2+), 3) HR-negative and HER2-positive (HR−/HER2+), and 4) HR-negative and HER2-negative (HR−/HER2−). Treatment modalities included surgery, radiation, chemotherapy (cyclophosphamide, epirubicin, and taxanes), trastuzumab, zoledronic acid, and hormone therapy (antiestrogen, tamoxifen, and aromatase inhibitors). Patients typically received a combination of these treatments. Noninvasive and early-stage invasive BCs (stages I and II) have more favorable prognoses than later-stage BCs (stages III and IV). Stage IV cancer spreads from the breast and nearby lymph nodes to other areas, resulting in a poor prognosis [10]. Treatment strategies were determined based on the patient’s molecular subtype and tumor stage and were associated with the patient’s prognostic factors. HR-negative patients did not receive hormone therapy, whereas all HER2-positive patients received trastuzumab according to the guidelines [11]. Because hormone therapy, trastuzumab, and chemotherapy may overlap with prognostic factors, these were not considered as independent factors in the analyses.

Patient population
We extracted information on BC-diagnosed women, aged 40 plus at time of diagnosis (earliest diagnosis year was 2006), from the nationwide Norwegian Cancer Registry. Patients were followed until the first CVD event (each of the six CVDs listed in the Introduction), death, or throughout 2020. Information on CVD events, identified according to International Classification of Diseases, 10th Revision codes, was obtained from the nationwide Norwegian Patient Registry [1213]. Data from both the Norwegian Cancer Registry and Norwegian Patient Registry were obtained through a common application process. Upon approval, the data were provided to the authors. The time of diagnosis was recorded in years, whereas the time of the CVD event was recorded in calendar days. For simplicity, the time of diagnosis was assumed to be the first day of the year. We included only patients who were available for CVD events in the calendar year following diagnosis. Therefore, the minimum follow-up time before a CVD event was 365 days.
We included patients with valid timelines, unilateral BC events in the diagnosis year, known tumor stages (I–III), surgery treatment, and radiation status, reducing the dataset from 45,334 to 36,778 patients (Figure 1 illustrates the data selection procedure). A small number of patients did not undergo any surgery. Most patients had stage I–III tumors. Patients with stage IV tumors were excluded because they generally have a poor prognosis and different survival patterns compared with those who have stage I–III disease.

Analytical approaches
Because some patients may die before experiencing a CVD event, death was treated as a competing risk. As our primary objective was to evaluate how CVD events are associated with prognostic factors (i.e., BC etiology), we used a cause-specific hazard modelling approach. However, the Fine-Gray alternative model has been described as a more appropriate alternative for making predictions in the presence of competing risks [14]. The Cox regression model for analyzing the time to an event assumed a proportional hazard (Cox-PH). However, this may not be the case, particularly in longitudinal data with extended follow-up periods [1516]. If underlying assumptions such as Cox-PH are not met, alternative non-parametric models such as RSFs and neural networks have become attractive alternative approaches. RSF is a well-established ensemble learning method for the survival analysis of time-to-event data that are possibly censored. It builds multiple decision trees, allows for nonlinear relationships, and does not assume proportional hazards [17]. Especially in cases where the factors have a complex non-linear relationship with the time-to-event outcome, the RSF becomes a good alternative to the Cox model. The RSF method is particularly feasible for determining and ranking the importance of each factor in a time-to-event outcome. Although our primary focus was on the impact of prognostic factors and their importance as obtained by the fitted RSF model, we also displayed the predicted cumulative incidence function (CIF) curves commonly generated by RSF models. These results complement the empirical cumulative incidence of HF and the competing risk of death. Additionally, we presented plots illustrating the nonlinear relationships between prognostic factors and time to HF based on the fitted Cox model.
We fitted a Cox model for the time to the first occurrence of each CVD, with death as a competing risk. Accordingly, patients who died without experiencing any CVD events were censored at the time of death. Associations between the molecular subtype groups and time to CVD, adjusting for age, tumor stage, and radiation therapy, were examined. All patients were assumed to have received chemotherapy, trastuzumab, zoledronic acid, and hormone treatment according to the guidelines [12]. Models were initially specified to include all factors, and an automatic model selection procedure based on the Akaike information criterion was applied to evaluate the optimal model fitting with respect to which factors to include, constraining the inclusion of the molecular subtype [1819]. Schoenfeld residuals were computed over time for each factor in the model. In a test for the independence between these residuals over time, small p-values indicate time dependence and, hence, a violation of the Cox-PH assumption [20].
We followed up with RSF analysis. The forest was based on constructing trees where for each tree we divided the dataset into a learning and test set of 67% and 33% of the data, respectively. From the test set, 1,000 bootstrap samples were drawn to construct risk trees. The tree was formed based on the log-rank split rule, where the algorithm randomly selected a subset of the factors at each node in the tree and split the node using the factor that maximized the splitting rule [17]. This was performed for 500 trees, and the average of this ensemble of trees was computed. The relative importance of molecular subtype, age, tumor stage, and radiation therapy was assessed, and CIFs were predicted over time for different groups according to molecular subtypes, age, and tumor stage given radiation therapy.
All analyses were performed using the R statistical program (R Foundation for Statistical Computing, Vienna, Austria) [21]. The survival and randomForestSRC packages were applied for Cox and RSF analyses, respectively [2223]. The Cox zph function in R was applied to test the Cox-PH assumption by computing the Schoenfeld residuals over time for each factor in the model [20].

Ethical statement, funding, and data access
This study was supported by The Dam Foundation (project number: 389462). All data were anonymized and made publicly available. This study was approved by the Regional Committees for Medical and Health Research Ethics (REK) (number 245944), which granted an exemption from the requirement for individual patient consent. Access to this dataset from Norwegian public registries can be obtained upon application and requires prior approval from the REK.

RESULTS

Patient population
A total of 7,039 patients (19.2%) experienced at least one CVD event during the study period, whereas 2,016 patients (5.5%) had at least one HF, and 3,401 patients (18.5%) died among those who did not experience at least one CVD.

Table 1 presents patient characteristics for all 36,605 patients, stratified by those who had no CVD and died (n = 3,401), experienced (n = 7,039) and did not experience (n = 26,165) at least one CVD, and experienced each of the six CVDs (n = 925 [MCI] to n = 2,982 [AF]). One patient might experience different CVDs. The median follow-up time until death, any first CVD, or the last year of observation (2020) was 6 years. The median follow-up time until death, any HF, or last year of observation (2020) was 7 years.
Most patients (63.0%) were HR+/HER2−, while only 6.8% and 3.6% were HR+/HER2+ and HR−/HER2+, respectively. The HR−/HER2− group accounted for 26.6% of patients. Compared with patients who did not experience any CVD, those who experienced at least one CVD were generally older and were less likely to have received radiation therapy (62.5% vs. 84.0%) and hormone therapy (59.8% vs. 75.3%) but more likely to have received chemotherapy (67.8% vs. 55.1%) and zoledronic acid (89.9% vs. 67.6%).
The incidence rates per 1,000 years for patients who experienced at least one CVD event and for patients who experienced each of the six CVDs are provided in Supplementary Table 1 in the supplemental material. The HER2-positive patients had the highest incidence of CVD per 1,000 years; 293 for HR+/HER2+ patients and 302 for HR−/HER2+ patients. The HER2-negative patients had lower CVD incidence: 262 for HR+/HER2− patients and 181 for HR−/HER2− patients.

Cumulative incidence of HF
The empirical cumulative incidences of HF and death over time, stratified by molecular subtypes, age groups, and tumor stage groups are shown in Figure 2. The CIF curves for HF in the HER2-positive patients crossed the CIF curves for mortality after approximately 13 years. The CIF curves for HF, especially death, increased with age and tumor stage.

Cox survival analyses
The optimal fitted Cox models for the time to the first AF, ES, and HF included the molecular subtype, age, tumor stage, and radiation therapy, while in the optimal models for the first occurrence of MCI, AP, and EPV, tumor stage factors were not a factor. In a test for the independence between the Schoenfeld residuals over time, small p-values indicate time dependence and, hence, a violation of the Cox-PH assumption, as shown in Figure 3, for the time to HF.

The MCI, ES, AF, AP, and EPC endpoints
In the model evaluating time to the first MCI, the Cox-PH assumption was satisfied for all molecular subtypes. No significant difference was observed between HR−/HER2− and HR+/HER2− patients (hazard ratio, 1.24; confidence interval (CI, 1.08−1.43). For the model for the time to the first ES, the Cox-PH assumption for the molecular subtype was also met; however, no significant differences were observed among the molecular subtype groups. For the time to the first AF, AP, and EPV model, the Cox-PH assumption was not satisfied for some molecular subtype groups.

The HF endpoint
Significant differences were observed between the molecular subtypes, as shown in Table 2. The Cox-PH assumption was not met for the molecular subtype factor, indicating that the association between molecular subtype and time to first HF varied over time. Additionally, this assumption was not met for the radiation therapy factor.
Compared with HR+/HER2− patients, HER+/HER2+ patients had a hazard ratio of 2.16 (CI, 1.85–2.52), whereas HR−/HER2+ and HR−/HER2− patients had hazard ratios of 2.44 (CI, 2.00–2.96) and 1.25 (CI, 1.13–1.38), respectively. The hazard ratio increased with age and tumor stage severity. Patients who received radiation therapy had a hazard ratio of 0.76 (CI, 0.68–0.84) compared with those who did not receive radiation. Table 2 presents the results from the cause-specific Cox regression model for time to death, treating time to HF as a competing risk.

Figure 3 displays the estimated HR over time. The blue line represents the HR-fitted smoothing spline that varies over time, with the gray area indicating the +/− 2-standard error band. The turquoise line shows the time-invariant hazard ratio estimated from the Cox model, assuming time-invariant coefficients, as shown in Table 2.
The hazard ratios estimated from the Cox model do not account for the time-varying associations between molecular subtype, radiation, and time to an HF event. An alternative to the Cox model is the accelerated failure time model, which does not assume a proportional hazard; however, it does also not adequately capture the temporal variation in the association between the molecular subtype and HF. More complex survival models that incorporate smooth splines, as shown in Figure 3, are of course possible, but one needs, as Nasejje and Mwambi points out, “to specify correct degrees of freedom, number and placement of the knot points, and order of the regression spline model (which could be quadratic, cubic, quartic, some combination of different orders, among others),” and they further point out that several models could fit the data just as well but estimate quite different hazards [24].

RSF analysis for time to HF
We performed an RSF analysis to estimate the cumulative cause-specific incidence (CIF) for time to HF and death for each tree and obtained the average over all trees. The importance of molecular subtype, age, tumor stage, and radiation therapy was ranked by computing the variable importance (VIMP), describing the impact of the factors on the prediction accuracy of the model.

Figure 4 displays the predicted CIF curves for HF and death for the four molecular subtype groups for patients in age groups 53–69 and 70–79 years old, with tumor stage I–III, who all received radiation therapy. For the time to HF, there was a difference in the CIF curves between HER2-positive and HER2-negative patients, where HER2-positive patients generally had CIF curves above the CIF curves of HER2-negative patients. Concerning time to death, we found the HR−/HER2− patients’ CIF curve to lie above the CIF curves for the other molecular subtype groups, except for HR−/HER2+ patients 53–69 years old with stage I tumors. Furthermore, the observed crossings of the CIF curves suggest a non-linear relationship between time to death and molecular subtype.

Table 3 presents HF and death-specific VIMP values (times 100). Factors with VIMP greater than 0.002 were considered important factors for survival prediction [25]. For time to HF, age was identified as the most predictive variable, followed by molecular subtype, radiation, and tumor stage. For time to death, age again showed the highest predictive variable, followed by tumor stage, molecular subtype, and radiation therapy.

DISCUSSION
The estimated CVD incidence rates per 1,000 years showed that patients in the HER2-negative group had a lower risk of CVD after BC treatment than those in the HER2-positive group, ranging from 181 for HR−/HER2− patients to 302 for HR−/HER2+). A previous study investigating cardiotoxic effects of chemotherapy drugs identified anthracyclines (1%–26% of patients), high-dose cyclophosphamide (7%–28%), and trastuzumab (2%–28%) as relevant factors, all commonly used in BC treatment [5]. Because all HER2-positive patients received trastuzumab and many also received chemotherapy, this combined regimen likely increased cardiotoxicity, resulting in more CVD events. These findings are consistent with those of Ngô et al. [26], who reported HER2-positive patients had a higher risk of cardiovascular mortality compared with HER2-negative patients.
For time to MCI Cox regression analysis, a significant difference was observed in hazard ratios between HR−/HER2− and HR+/HER2− patients (hazard ratio, 1.24; CI, 1.08–1.43), indicating that the risk of developing a first MCI was 24% higher in HR−/HER2− patients than in HR+/HER2− patients. Tumor stage was not a significant factor in this model. No significant differences were observed in the risk of ES among different molecular subtypes. For the other endpoints, risk differences were observed among some molecular subtype groups; however, the Cox-PH assumption was not satisfied. Finally, for the HF endpoint, significant differences were observed among molecular subtype groups, but the Cox-PH assumption was not satisfied.
The association between the time to first occurrence of HF and molecular subtype clearly varied over time, as illustrated in Figure 3. The Cox-PH assumptions were not met, and the estimated HRs presented in Table 2 should be interpreted as average-over-time hazards. Based on these results, HER2-positive patients had a higher risk of developing HF than HER2-negative patients. This pattern was also observed in the predicted CIF curves from the RSF analysis (Figure 4), in which the HER2-negative and HER2-positive patients were separated into two groups. Both the Cox regression and random forest model analyses revealed a clear distinction between HER2-positive and HER2-negative patients with respect to the risk of developing HF following BC diagnosis and treatment. Additionally, HR−/HER2− patients showed a slightly increased risk of HF compared with HR+/HER2− patients. Because HR−/HER2− patients do not receive trastuzumab, they are likely to undergo intensive chemotherapeutic treatments (including cyclophosphamide, epirubicin, and taxanes), accounting for this increased risk [7]. Thus, different molecular subtype groups lead to differences in treatment. This results in differences in risk of HF, and these differences change over time.
The time-varying risk patterns of CVD between HER2-positive and HER2-negative BC patients can be explained by distinct treatment-related mechanisms, biological factors, and patient characteristics. For HER2-positive patients, targeted therapies, such as trastuzumab, pertuzumab, and trastuzumab emtansine, are associated with early-onset cardiotoxicity, typically within the first year of treatment. HER2-positive patients are typically closely monitored with frequent cardiac surveillance, allowing for the early detection of cardiac issues [27]. HER2-negative patients typically tend to be older and have higher baseline CVD risk in terms of, for example, hypertension and diabetes. HER2-negative patients might not receive routine cardiac surveillance unless symptoms are expressed, resulting in delayed recognition of cardiovascular issues.
An advantage of using the RSF approach is that it provides a ranking of factors, molecular subtype, age, tumor stage, and radiation therapy, with respect to their predictive importance for the time to first occurrence of HF. Age was identified as the strongest and most important factor, which was not unexpected [28]. The molecular subtype group was identified as the second most important variable, reflecting the treatment decisions associated with cytotoxicity in these patients.
Because this was a population-based analysis, the study was not affected by observational bias. The dataset included information on women diagnosed with BC in Norway during the study period.
Because the Cancer Registry data contained only annual dates, whereas the Patient Registry data on CVD events included exact dates, certain assumptions regarding the date of BC onset were necessary. Furthermore, information regarding the dose, frequency, or time of treatment was unavailable. However, it is reasonable to assume that the treatments were administered in accordance with international guidelines. The Norwegian Cancer Registry, operated by the Norwegian Institute of Public Health, is mandatory and considered to have very high completeness and validity, with multiple independent data sources contributing to its robustness. Therefore, the data was considered reliable and of good quality.
In conclusion, the relationship between the time to first CVD and molecular subtype differs depending on the specific CVDs, highlighting the need to evaluate each CVD endpoint independently. Our analysis revealed that the prognostic significance of factors such as molecular subtype can change over time, which has not been well described previously. RSF analyses can capture complex relationships between predictors and survival without requiring the Cox-PH assumption and can provide important insights into how patient and treatment factors are associated with the time of CVD events.