Separate and combined analysis of successive dependent outcomes after breast-conservation surgery: recurrence, metastases, second cancer and death

Several randomized trials have shown that breast-conserving surgery (BCS) plus locoregional radiotherapy is an appropriate method of primary therapy for the majority of stage I and II breast cancers and may be preferable to radical mastectomy because it provides survival equivalence while conserving the breast [1]. The main concern for both physicians and patients is therefore the risk of local or metastatic recurrence in the conserved breast or in the contralateral breast.

The risk of recurrence (local, locoregional or contralateral), distant metastases (DM) or second primary cancer after breast cancer and their predicting factors has already been documented in the literature [2–6]. In most applications, Cox proportional hazard models are fitted and prognostic factors are studied for each specific event [3]. The most straightforward approach in these recurrent event settings is simply to count only the first occurrence of an outcome in a subject [5]. This analytical strategy is straightforward, and it avoids the methodological complications that can occur if a first event affects either the risk of subsequent events or compliance with ongoing treatment [7]. However, consideration of only first events is not satisfactory to evaluate the natural history of a disease and the effects of therapeutic interventions, furthermore, this is inefficient because it does not utilize all the available information. The potential benefits in terms of events prevented by a treatment can be underestimated if only the first failure is considered.

The solution could be to include as time-dependent covariates the different intermediate events which may affect the patient's event of interest and to study their prognostic impact [4, 8]. Yet, interpretation of these time-dependent covariates can be difficult in practice. In this setting of semi-competing risk events, Cox proportional hazard models are simply used for each of the events [9], yet they are unable to estimate the association between these multiple events.

This article hypothesizes that Cox proportional hazards models are not the most appropriate for studying the evolution of breast cancer-related events occurring in a patient. For a more detailed interpretation, one should consider all events. The latter are dependent, for example ipsilateral breast tumor recurrence may be associated with subsequent DM and worse survival [4, 10]. A frailty model is proposed to identify the prognostic factors of breast cancer recurrences, metastases and second primary cancer among women treated by BCS for a primary invasive breast cancer. Multiple events in the same subject are considered so it is important to control for the correlation of recurrent events of different types. The frailty model is suitable for studying these recurrent events and prevents biased inferences [11]. In cancer studies, population heterogeneity can be particularly important because of different exposures to carcinogens, different genetic predispositions [12] and the considerably varying speed with which the disease evolves. This population variability will be expressed by some women remaining free of recurrences throughout follow-up and by others having frequent events. Such frailty models could identify the risk of metastases and the need for an adjuvant systemic treatment after a local recurrence. By examining the association between various breast cancer-related events, it may be possible to develop preventive interventions. Furthermore, the statistical power is increased when all available outcomes are considered together instead of using separate statistical analyses for each subsequent outcome.

The number of events ranged from 0 to 4 (average = 0.62) per patient. The first failure type was mostly LRR (n = 134, 14.6%). Second failures were most commonly death from breast cancer (n = 76, 22.4%) for all first failure types. The most frequent event observed over the whole follow-up was metastases, which was observed in 188 of the women (13.9% of the total number of events observed).

Analysis of each event separately

The results of the five separate analyses for each event of interest are shown in Table 3.

Characteristic	Locoregional recurrences		Contralateral breast cancer		Distant metastases		2nd primary cancer other than breast		Breast cancer death
	RR (95% CI)	p	RR (95% CI)	P	RR (95% CI)	P	RR (95% CI)	p	RR (95% CI)	p
Histological size of the tumor (for 1 mm increase)	1.02 (0.99-1.06)	0.17	1.02 (0.97-1.06)	0.54	1.02 (1.00-1.05)	0.08	1.01 (0.93-1.10)	0.75	1.02 (0.99-1.05)	0.27
Extensive intraductal component (EIC)
- No (ref)	1		1		1		--		1
- Yes	3.34 (0.89-12.5)	0.07	0.99 (0.13-7.46)	1.00	1.55 (0.34-7.06)	0.57	--	--	1.43 (0.20-10.8)	0.73
Grade
- Grade1(ref)	1		1		1		1		1
- Grade 2	1.84 (1.11-3.07)	0.02	1.45 (0.78-2.70)	0.24	2.68 (1.58-4.53)	< 0.001	0.68 (0.23-2.01)	0.48	3.64 (1.74-7.60)	< 0.001
- Grade 3	2.12 (1.17-3.84)	0.01	0.94 (0.40-2.21)	0.88	3.39 (1.91-6.02)	< 0.001	1.32 (0.37-4.75)	0.67	5.54 (2.58-11.9)	< 0.001
# positive nodes
- 0 (ref)	1		1		1		--		1
- 1-3	2.60 (1.30-5.19)	0.007	1.17 (0.34-3.97)	0.81	2.41 (1.29-4.47)	0.006	--	--	2.92 (1.50-5.68)	0.002
- ≥4	4.40 (1.80-10.7)	0.001	0.66 (0.13-3.45)	0.62	6.91 (3.41-14.0)	< 0.001	--	--	8.03 (3.71-17.4)	< 0.001
Lympho-vascular invasion
- No (ref)	1		1		1		1		1
- Yes	1.29 (0.84-1.98)	0.24	0.41 (0.19-0.90)	0.03	1.47 (1.04-2.08)	0.03	1.40 (0.51-3.83)	0.51	1.27 (0.87-1.91)	0.26
Surgical excision margin status
- Clear (ref)	1		1		1		1		1
- Invasion	0.95 (0.64-1.43)	0.81	2.36 (1.40-3.98)	0.001	1.05 (0.74-1.48)	0.78	0.57 (0.19-1.73)	0.32	1.12 (0.75-1.67)	0.59
Progesterone receptor status
-Positive(ref)	1		1		1		1		1
- Negative	1.26 (0.84-1.88)	0.26	1.21 (0.70-2.09)	0.49	1.32 (0.93-1.86)	0.12	0.33 (0.11-1.03)	0.06	1.90 (1.26-2.87)	0.002
Treatment combination
- Radiotherapy (ref)	1		1		1		--		1
- Radio- and chemotherapy	0.23 (0.11-0.50)	< 0.001	1.22 (0.33-4.55)	0.77	0.52 (0.27-1.02)	0.06	--	--	0.56 (0.27-1.16)	0.12
- Radio- and hormonal therapy	0.35 (0.16-0.77)	0.009	0.78 (0.21-2.91)	0.71	0.59 (0.30-1.13)	0.11	--	--	0.48 (0.24-0.99)	0.05
- Radio- and chemo- and hormonal therapy	0.06 (0.01-0.49)	0.008	0.74 (0.07-7.60)	0.80	0.33 (0.11-0.99)	0.05	--	--	0.28 (0.07-1.06)	0.06

-- variables not included in the model

They suggest that the risk of LRR was higher for women with higher tumor grade, with an extensive intraductal component or with nodal invasion. In patients treated with a combination of radiotherapy, chemotherapy and hormonotherapy, less locoregional occurrences were observed.

The risk of contralateral breast cancer was increased in patients showing involvement of excision margins after BCS (RR = 2.36, 95%CI = 1.40-3.98). By contrast, lympho-vascular invasion was inversely associated with contralateral breast cancer.

With regard to factors influencing distant metastases, it was found that grade, nodal invasion and lymphovascular invasion were highly prognostic factors of distant metastases. In patients treated with a combination of radiotherapy, chemotherapy and hormonotherapy less metastases were observed.

No predictive factors except the progesterone receptor status were associated with the risk of developing a new second primary cancer other than breast carcinoma, but a poor statistical power may explain these results (only 30 new second primary cancers were observed).

Grade and number of positive nodes were also associated with breast cancer death. Progesterone receptor positivity decreased the risk of death from breast cancer.

The following factors were not associated with any of the five events of interest and were not retained in analyses in Table 3 and 4: histological type, family history of cancer or breast disease, parity, estrogen receptor status.

	Model 1	Model 2	Model 3
Characteristic*	RR (95% CI)	RR (95% CI)	RR (95% CI)
Histological size of the tumor (for 1 mm increase)	1.03 (1.01-1.06)	1.02 (1.00-1.04)	1.02 (1.00-1.04)
Extensive intraductal component (EIC)
- No (ref)	1	1	1
- Yes	1.91 (0.58-6.24)	1.44 (0.58-3.59)	1.21 (0.49-3.01)
Grade
- Grade 1 (ref)	1	1	1
- Grade 2	2.06 (1.44-2.95)	1.77 (1.30-2.40)	1.69 (1.24-2.29)
- Grade 3	3.19 (2.08-4.87)	2.41 (1.69-3.43)	2.32 (1.64-3.29)
# positive nodes
- 0 (ref)	1	1	1
- 1-3	2.82 (1.62-4.91)	2.33 (1.52-3.57)	2.33 (1.52-3.55)
- ≥4	6.16 (3.15-12.04)	4.60 (2.74-7.72)	3.97 (2.36-6.66)
Lympho-vascular invasion
- No (ref)	1	1	1
- Yes	1.10 (0.80-1.52)	1.08 (0.85-1.38)	1.02 (0.80-1.30)
Surgical excision margin status
- Clear (ref)	1	1	1
- Invasion	1.36 (1.02-1.82)	1.21 (0.96-1.52)	1.23 (0.97-1.54)
Progesterone receptor status
- Yes (ref)	1	1	1
- No	1.45 (1.09-1.93)	1.43 (1.14-1.80)	1.46 (1.16-1.83)
Treatment combination
- Radiotherapy (ref)	1	1	1
- Radio- and chemotherapy	0.50 (0.27-0.91)	0.55 (0.35-0.87)	0.59 (0.37-0.93)
- Radio- and hormonal therapy	0.36 (0.19-0.66)	0.50 (0.31-0.81)	0.54 (0.34-0.86)
- Radio- and chemo- and hormonal therapy	0.22 (0.08-0.63)	0.29 (0.12-0.67)	0.31 (0.13-0.72)
Number of prior events (as a quantitative variable)	-	1.97 (1.62-2.41)	0.96 (0.67-1.41)
Type of prior events:
- no previous event (ref)	-	-	1
- locoregional	-	-	4.25 (2.51-7.21)
- contralateral breast cancer	-	-	1.44 (0.75-2.74)
- distant metastases	-	-	3.94 (2.23-6.96)
- 2nd primary other than breast			- (-)
Variance of the random effect	1.42 (0.24)	0.38 (0.20)	0.35 (0.17)
Akaike Information Criterium (AIC)	1518.98	1494.10	1461.69

* Analyses parametrically adjusted for the characteristics mentioned and non-parametrically adjusted for age

Figures 1 and 2 show the baseline recurrence hazards generated from adjusted models, according to the rank of the recurrence or the type of the recurrence. Figure 1a shows that the first two baseline age-specific recurrence hazards are higher for young women (less than 40 years old). It also shows that at any age the risk of failure (recurrence, contralateral breast cancer, metastases, 2nd primary other than breast or death) increased successively with the number of prior events. For instance, the risk of developing a third event if a patient had already developed a second event was higher than the risk of developing a second event if she had already developed a first event. Figure 1b illustrates the four separate time-specific recurrence hazards and shows that the risk of developing a first relapse event is fairly stable across time since surgery. The time to first event is longer than the other gap times (i.e. time between two successive events) and can therefore occur at an older age.

Figures 2a and 2b show the 5 type-specific hazards (for recurrence, contralateral breast cancer, metastases, 2^nd cancer other than breast or death from breast cancer) without taking into account the rank of failure. The different types of events follow different patterns across the time intervals. The risks of developing LRR or DM is higher than the risks of developing other events at any time and at any age. A greater risk of metastases was observed during the first 10 years after surgery (see Figure 2b) but the other types of failure occurred at a constant rate over at least 10 years. After that time, the risk of failure decreases progressively.

This study uses a novel approach of analyzing together the different rates of failure following successive events for women with breast cancer treated with BCS. The risk of relapse of any type was higher for young women (less than 40 years old) and is higher in the first 10 years after surgery. Traditionally guidelines for follow-up of breast cancer patients concentrate on the first 3-5 years, with either reduced frequency of visits or discharge after this. Montgomery et al. [18] show that the basic assumptions behind these generally-accepted guidelines for follow-up may be incorrect. They suggest that if the central aim is the detection of treatable relapse, as stated in the guidelines, then there is no justification for focusing on the first 2-3 years after initial therapy, as treatable relapses occur at a constant rate over at least 10 years. Our results are in accordance with these results and the implications for cancer care should be carefully considered.

As this was a prospective observational study in a single center, there are potential limitations to our results. Principal among these limitations could be a bias for patient selection or for their follow-up. In terms of patient selection, we noted a median age (57 years) lower than the median age incidence in France [19] (61 years), but we selected only invasive breast cancer. The potential bias for follow up in these types of studies appears not to be the case in this study as the cancer center employs a rigourous clinical follow-up procedure and we administered a survey to physicians resulting in less than 2% losses to follow-up [13].

This frailty model allows the analyst to estimate simultaneously and distinctly the effects of both unobserved heterogeneity through a subject-specific random effect and event dependence through a time-dependent covariate. It also helps us to diagnose the source of correlation in the data by comparing the magnitude of the estimated effects across the models, i.e. the adjustment to some covariates will completely explain the unobserved heterogeneity. In this typical repeated events context, the usual Cox model was both biased and inefficient and the five separate analyses were uninformative with regard to heterogeneity across individuals and event dependence. Furthermore, the combined analysis increased the statistical power to such an extent that the associations of some prognostic factors with any event were revealed (ex: histological tumor size in Table 4). No separate analysis was able to do this. Since the risk of any recurrence was higher for young women and in the first 10 years after surgery, this justifies the need for increased surveillance of women in the first years following the diagnosis of cancer. Tai et al. [20] recently proposed a marginal approach in order to utilize fully all the available information on event times and not only the information on the first event that occurs. But, their approach does not measure the dependence among the multiple events and the influence of prior events on future failures is not their focus. Recently using a multivariate multistate model, de Bock et al. [21] found that patients with a LRR have a more than threefold increased risk of developing DM as compared to patients who develop no LRR. As with the frailty model, the strongest point of their approach is that all the data are summarized in one model instead of presenting many separate analyses. Presenting the many separate analyses will lower the power of the estimated effects or may result in false positive findings. However such a multistate analysis can only be performed on a large cohort of patients with a long follow-up time (with enough metastatic events after the occurrence of LRR). The frailty model is less restrictive on the number of events.

At the end of the study, 209 patients had died, among them 133 had died from their breast cancer. Deaths could be related to their health status and more specifically to breast cancer recurrences (local or distant). Thus, breast cancer death is likely to be an informative type of censoring, which means that those individuals who are censored by death are not as likely to have the subsequent event of interest as those who remained in the study. In order to examine this hypothesis we used a joint frailty model (using the R package "frailtypack") to analyze recurrent observations of breast cancer with death from breast cancer [22, 23]. For instance, a patient recurrence rate may be positively correlated with its death rate. The effects of prognostic factors were unchanged in this joint frailty model. Therefore, censoring by death from breast cancer appears to be non-informative for recurrence after adjustment. This suggests that selection by death from breast cancer in this study did not importantly bias our findings.

In this study we proposed to control for the type of event by adjusting for this variable in the frailty models (see Table 4, model 3). We then observed that patients with a LRR had a more than fourfold increased risk of developing a new event as compared to patients who did not yet develop any event. However in this approach the prognostic factors had the same effects for the different types of event. Another approach could be to use a joint frailty model with one hazard function for each type of event [24]. This would allow us to evaluate the prognostic impact for each type of event and so for each clinical etiology. However we have not considered this extension in the present application, due to the increased computational complexity of the problem.

Using breast cancer data from a population-based cancer registry, Cheung et al. [25] illustrated the features of Cox models using two time scales (time since diagnosis and age). Using time since diagnosis as the time scale, a younger age at diagnosis was associated with a lower mortality, while using age as the time scale gave the opposite result. Age was chosen as the basic time scale in our analysis, so no parametric assumptions were made in the analyses and the hazard function of the time of onset of any recurrence was the age-specific incidence of the recurrence. We also analyzed parametrically the effect of age by introducing age as a covariate in a proportional hazard model (i.e. with a log-linear relationship). The results were quite similar in models adjusted for age with four classes coded by three binary variables (40-49, 50-69, 70 and over versus 29-39 years). Age was associated with the risk of locoregional recurrence or metastases with a higher risk for younger patients. However, if the true effect of age is far from this parametric log-linear assumption, this may lead to spurious results.

Our study seems to show a previously unreported unexplained protective effect of lympho-vascular invasion on the risk of contralateral breast cancer. The epidemiology of contralateral breast cancer is complex to study for a number of reasons. Firstly, the patient has received treatment that may modify the risk of second cancer, and she is under intense medical scrutiny and self-observation [26, 27]. Secondly, the treatment and follow-up procedures for breast cancer have changed extensively over the years and several definitions of contralateral breast cancer were observed. This may contribute to the incoherent results reported in the epidemiology of contralateral disease [27, 28]. Comparison of genomic profiles would allow differentiation between metastatic lesions and new primary tumor, but this information was not available to us. Microarray studies have shown that distinct molecular profiles can be useful for classifying tumors and predicting outcome [29].

Tailored systemic treatment making use of ER and PR receptors and HER2 is now current practice [30]. However these prognostic factors were missing in our study. Only more recent studies will include the newer predictive and prognostic factors, but these studies will not allow us to perform such long-term analyses due to their shorter follow-up.

In conclusion, the frailty model approach allows the simultaneous analysis of all successive events intervening after a first breast cancer. This methodology in conjunction with conventional methods represents a valuable tool for describing the natural history of breast cancer and the possible effects of treatment. This methodology issue may also be relevant to studies of other tumor types with recurrences. This approach can also, if data require it, incorporate the effects of interventions which are performed after each event. Response to treatment after each relapse may be an important factor to predict new relapses; this is an important issue which is not considered in the models for recurrent event data.

We consider models in which the hazard function partly depends on an unobservable random variable thought to act multiplicatively on the hazard, so that a large value of the variable increases the hazard. These models, called frailty models are an extension of the classical Cox proportional hazard model [15].

where λ ₀(t) is the baseline hazard function; X _ij = (X _1ij ,..., X _pij )' denotes the covariate vector for the j^th recurrence of individual i, β is the corresponding vector of regression parameters, and the Z _i 's are unobserved random variables (the frailties). It is assumed for mathematical convenience that the Z _i 's are independently and identically distributed from a gamma distribution with mean 1 and unknown variance θ at time of origin.

Different timescales can be used [31]. The timescale that is most often used is the gap time: after an event, the subject starts again at time 0 and the time to the next event corresponds to the number of days that it takes to experience the next event. An alternative timescale is the calendar time, also called the counting process approach [32] which keeps track of time since randomization. The duration of the time at risk for an event corresponds to the duration of the time at risk in the gap time representation. However, the start of the at-risk period is not reset to 0. A subject is not considered to be at risk for the k^th event until after the (k-1)^th event. A particular subject has different periods at risk during the total observation time. If there are n _i at-risk periods for patient i, then the complete information for patient i can be represented by n _i triplets $(Y_{i11},Y_{i12},{\Delta }_{i1}),\mathrm{...},(Y_{i{n}_{i}1},Y_{i{n}_{i}2},{\Delta }_{i{n}_i})$ where, for the j^th triplet, Y _ij1 is the start of the j^th at-risk period, Y _ij2 is the end of the j^th at-risk period, Δ _ij is the censoring indicator and Y _i11 = 0. If age is chosen as the basic timescale, the hazard functions estimated are for instance the age-specific recurrence incidences or the age-specific mortality. A subject is considered at risk from her age at surgery (= Y _i11) to her age at censorship or failure.

The expression of the full marginal likelihood associated with this frailty model is used. In most situations it is reasonable to expect smooth baseline hazard functions, piecewise constant modeling for the hazard functions often being unrealistic. To introduce such a priori knowledge, we penalize the likelihood. The expression of the full marginal likelihood is developed in expression 4 of Rondeau et al. [23].