A multi-state model incorporating estimation of excess hazards and multiple time scales

Marked improvement in cancer patient survival in recent decades has increased the importance of studying treatment side-effects[1], i.e., late effects. For example, chemotherapy and radiotherapy increase the risk of both secondary malignancies and cardiovascular disease (CVD)[2, 3]. As therapies are developed and adapted to maximise the chances of cure while reducing risk of treatment toxicity, it is of interest to simultaneously estimate real-world morbidity and mortality. Today, a major challenge in the field of late effects lies in estimating the proportion of the observed morbidity and mortality that can be ascribed to the underlying cancer and its treatment. Hence, the standard approach is to simply estimate the cumulative observed incidence and mortality[4]. This article proposes a population-based approach to this problem bringing together statistical methods for relative survival and multi-state modelling. The former allows estimation of excess incidence, which under certain assumptions corresponds to the incidence above and beyond that expected in the absence of cancer and treatment thereof, and the latter enables investigation of patient trajectories going from an initial state (diagnosis), potentially via intermediate states (e.g., non-fatal illnesses such as CVD), to final absorbing states (different causes of death). Our focus is on the illness-death model where the illness state can be experienced also in the absence of cancer (unlike recurrence, which only occur in cancer patients). In addition, we know only when patients enter each state and not whether it was due to the cancer. We are interested in partitioning the probability of being in the illness state into that due to the cancer and its treatment, and that which would be expected in the absence of cancer.

Previous attempts to combine relative survival with multi-state modelling has proposed a Markov piecewise constant rate model incorporating estimation of excess mortality[5, 6]. This method has been further generalised by Gillaizeau et al. to include semi-Markov models[7]. In these papers, the relative survival methods are applied to absorbing states alone, i.e., excess mortality can be estimated but not excess incidence. Moreover, focus has been on relative effects rather than predictions. Our suggested approach enables partitioning of both excess morbidity and mortality in a Markov or semi-Markov setting, and we are unaware of existing methods for doing this.

In population-based cancer patient survival, relative survival is a commonly used summary measure of patient survival as it captures both the direct and indirect contribution of a cancer diagnosis on mortality by comparing the observed survival of the patients to the expected survival in a comparable cancer-free population[8, 9, 10]. In the standard setting, the outcome is all-cause death and the mortality rates expected in the absence of cancer are typically retrieved from publicly available population mortality tables. However, the same approach can be used to estimate excess incidence rates[11]. As such, the excess incidence of some disease of interest can be interpreted as the portion of incidence attributable to the underlying cancer and its treatment, given that the expected rates are appropriate for the population diagnosed with cancer.

Lambert et al. have previously developed methods for estimating cumulative probabilities of death parametrically in a relative survival framework [12]. Eloranta et al. extended these methods to enable partitioning of the excess mortality associated with cancer into component parts such as excess cardiovascular mortality and remaining excess cancer mortality[13]. This has been applied to study temporal trends in excess mortality due to diseases of the circulatory system (DCS) among Hodgkin lymphoma (HL) patients[14]. Further extensions of these methods have been applied to study excess DCS incidence in HL patients[11]. In the latter paper, the expected incidence rate was modelled using a constant rate model with two time scales (calendar year and age) with twoway interactions. In the mentioned papers, the cumulative probabilities were calculated using numerical integration, in the presence of competing risks. However, considering morbidity as an absorbing state, as in the application on HL and excess DCS incidence, prevents identification of fatal complications, which in turn makes comparisons with previous studies on excess mortality more difficult.

Multi-state models are a generalisation of competing risks models and one way to study patient trajectories[15, 16]. The most commonly used method to estimate transition rates in a multi-state model is using a semi-parametric approach, such as the Cox model. In a recently published study, Crowther and Lambert[17] suggested a framework for transition-specific parametric distributions in multi-state models. In general, parametric survival models have several advantages over semi-parametric models, including the potential of modelling more than one timescale, more efficient incorporation of time-dependent effects, ease of prediction (both in- and out-of-sample), and greater flexibility in quantifying absolute risks. In addition, estimation of uncertainty is easier, and they enable extrapolation of survival over the full life span (although this should always be done with caution). In the previously proposed method[17], a separate parametric model can be fitted to each transition, which allows the use of transition-specific distributions. The transition probabilities are calculated using a simulation approach[18, 19]. Due to the flexibility of the simulation approach, other quantities of interest can be easily predicted, such as the average time spent in each state.

In survival analysis, there may be more than one plausible timescale on which to quantify an event rate[20]. An obvious example is how incidence, prevalence, and mortality due to various cancers vary as functions of attained age and calendar time. Acknowledging this may introduce transition-specific and/or multiple timescales, which are rarely investigated since the computational complexity increases substantially. Work in this area has been sparse, particularly moving to competing risks and the more general multi-state setting. Lee and Fine[21] developed cause-specific Cox models in a competing risks setting, allowing cause-specific timescales, where attained age was the underlying timescale for one cause of death and time since diagnosis was assumed for the competing cause of death. Wolkewitz et al.[22] developed a two timescale Cox model applied to modelling risk of infection among patients admitted to Intensive Care Units (ICU), with discharge and death as competing risks, using time since entry to the ICU and calendar time as timescales. In a parametric framework, the most substantial contribution to the field is by Iacobelli and Carstensen[16], who provided a general approach for estimating multi-state models with multiple timescales, centred on time-splitting and estimating rates within a parametric spline-based modelling framework using Poisson regression. Predictions of transition probabilities were obtained using simulation, and confidence intervals calculated using a parametric bootstrap.

Motivated by the need to simultaneously investigate excess morbidity and mortality among cancer survivors, this study proposes one solution to several methodological challenges related to the clinical application. Starting from an illness-death model (with cancer diagnosis as the initial starting state) we incorporate estimation of excess hazard rates into the transition model for the illness state. Moreover, in the expected hazard rate model we allow for multiple time scales (calendar time and age) while remaining transitions have one single time scale (time since study entry or time since entry to current state). Generalisation beyond the illness-death model is easily achieved. The approach is illustrated on a study of excess incidence and mortality from diseases of the circulatory system (DCS) among patients treated for HL.

The article is structured as follows. In Section 2 we give a brief overview of recently developed methods for multi-state modelling in the cause-specific framework, after which we describe our proposed approach for each transition. Section 3 provides an illustrative example on patients diagnosed with HL and excess DCS morbidity and mortality. We conclude with a discussion and suggest possible extensions in Section 4.

We begin by conceptually describing the method proposed by Crowther and Lambert[17]. We define a global vector of available covariates, $\bm{X}$, since the same covariates may be used in multiple models. The fitted transition-specific models can be as complex as required (with respect to functional form of covariates, interactions, and time-varying effects) and one can use different parametric models for different transitions (e.g., a flexible parametric for one transition and a Weibull model for another). Once the set of transition models has been chosen, transition probabilities are calculated for a covariate pattern of interest. The details have been extensively described elsewhere[17] but in short, point estimates are obtained by simulating a large sample of patient trajectories ($n$) through the fitted multi-state model. Confidence intervals are calculated from repeated simulations ($m$), each drawing from a multivariate normal distribution with the estimated parameter vector as mean and the variance-covariance matrix (of the parameters) as variance, i.e., a parametric bootstrap.

Figure 1 shows a representation of our suggested extension. While a standard illness-death model consists of three states (healthy, ill, dead), we partition the illness state into two unobserved states - excess and expected illness (dashed lines). We further have two death states, death following illness, and death without having experienced the illness, which allows us to quantify the impact of the illness under study on mortality. In total, we have four states, of which one is partitioned into two components. Patients start out in state 1; alive after a diagnosis of cancer. Patients might then experience some illness (state 2), which can be either related to the cancer and its treatment (“excess”) or not (“expected”). The transition rate to the expected illness component of state 2, $h_{1}^{*}(t+a,t+c)$, is modelled using population records of disease incidence with attained age ($t+a$) and calendar time ($t+c$) as time scales. The transition to the excess illness component of state 2, $\lambda_{2}(t)$, is estimated using a parametric excess hazards model, incorporating the modelled expected illness rate. Patients might die prior to experiencing the illness under study (state 3). This mortality rate, $h_{3}(t)$, can be modelled as all-cause or partitioned into cause-specific mortality models. The time scale for both the excess incidence model and this mortality model is time since cancer ($t$). Some of the patients experiencing the illness (excess or expected) will transition to the absorbing state of death (state 4). This transition, $h_{4}(t-u)$, can be modelled as an all-cause or cause-specific mortality model, with either one time scale (here we illustrate using time since illness, $t-u$) or multiple time scales (a natural extension would be to add $t$). Any of the models used for transitions into states 3 and 4 can be extended to excess hazards models, given the appropriate expected rate information.

We illustrate the suggested approach with an application to the modelling of excess incidence and mortality from diseases of the circulatory system (DCS) among patients treated for HL. HL is often pointed out as a model disease in terms of improvements in survival, going from a fatal disease in the mid-1900 to around 90% of young HL patients today being cured from their lymphoma[26, 27]. These improvements are primarily due to major advances in the management of these patients, including more accurate radiotherapy, effective multi-agent chemotherapy, and improved staging procedures. With more patients surviving their HL, research focus has somewhat shifted towards understanding and reducing treatment toxicity. The two major complications with potential fatal outcome are secondary malignancies and DCS.

From the Swedish Cancer Register, we identified 4,479 patients diagnosed with HL between 1985 and 2013, aged 18-80 at diagnosis. This cohort has been described in detail previously[11]. Patients were followed until death, emigration, or end of study period ($31^{\text{st}}$ of December, 2014), whichever came first. Incident cases of DCS were linked to the HL patients using the Swedish Inpatient Register. Only events with DCS as the main diagnosis were considered. At the end of each individual’s follow-up, 2,556 patients had not experienced any of the secondary events, 978 patients died due to any cause prior to experiencing DCS, and 945 patients had been diagnosed with DCS. Among these, 425 were alive and living with DCS at the end of study, and 520 patients had died (from any cause of death) (Figure 2).

For DCS, no population incidence table with sufficiently detailed data is publicly available in Sweden. Instead, a cohort of approximately 10 million individuals (Swedish residents who were either born after 1931 and residing in Sweden after 1961, and their parents, or who were included in a census 1960–1990) was constructed using several population-based registers, as a representation of the Swedish population. Data were stratified by calendar year, age, and sex, and the number of DCS events and amount of person-years at risk were calculated. The expected DCS rate was modelled using Equation 2, with 5 degrees of freedom for each time scale and incorporating sex.

For modelling all remaining transitions, flexible parametric survival models were used with year of diagnosis, age at diagnosis, and sex included as baseline covariates. Year and age were parameterised using restricted cubic splines with 5 degrees of freedom in the models for excess DCS incidence and mortality before DCS. For post-DCS mortality, 3 degrees of freedom were used due sparseness of data. No parameter effects were shared between models. For simplicity, only main effects of covariates were considered, however, interaction terms could easily be incorporated.

The transition from HL to excess DCS was modelled with time since HL diagnosis as the time scale. The predicted expected rate was incorporated into the likelihood to enable estimation of the excess rate. For the transition from HL to dead before DCS, time since HL diagnosis was used as time scale and all-cause mortality was the outcome of interest. Patients with DCS identified from their cause of death certificate ($n=100$) were given that date as their diagnosis date and the following day as their death date. Because of this, the all-cause mortality before DCS implicitly did not include any deaths caused by DCS. Death following a DCS diagnosis was modelled with time since DCS as the underlying time scale (clock-reset).

Figure 3 shows the transition rates for male patients aged 30 and 60 year old at diagnosis. The expected DCS rate is presented as a function of calendar time, while the remaining rates are for patients diagnosed in 1995, and plotted over time since diagnosis (cancer or DCS). The expected rates declined slightly following the mid-1980s, and the rate among younger patients was very low throughout. The excess DCS rate and both mortality rates were all high early during follow-up.

Using the transition-specific model estimates, point estimate of the transition probabilities were calculated by simulating 1,000,000 subjects through the multi-state model. Figure 4 shows the stacked transition probabilities over time since HL diagnosis for male patients diagnosed in 1995 at ages 30 (left panel) and 60 (right panel) years at diagnosis. For 30-year-olds, the probability of remaining in the HL state at 15 years after cancer diagnosis was approximately 80% whereas for 60-year-olds this probability was just above 20%. The probability of being diagnosed with DCS within the first 15 years as expected in the absence of HL and its treatment, and still being alive, was very small for younger patients. The probability of death (before or after a DCS diagnosis) was expectedly much larger for older patients.

Although stacked graphs provide a useful tool for illustrating in which state the patients are at different points in time, they lack information on uncertainty. As no parameters were shared between transition models, re-sampling to obtain confidence intervals was done using the transition-specific coefficient vectors and associated variance-covariance matrices, with $n=10,000$ and $m=1,000$. Figure 5 shows the transition probabilities with associated confidence intervals for HL patients diagnosed in 1995 at ages 30 and 60 years.

The fitted models can further be used to predict other clinically relevant measures, such as differences in transition probabilities across varying covariate patterns. Figure 6 illustrates the differences in transition probabilities between females and males diagnosed in 1995, aged 30 (top panel) and 60 (bottom panel) years.

Another measure of clinical importance is the proportion of excess DCS (calculated as the probability of excess DCS over the probability of any DCS), which allows us to identify trends in the contribution of excess DCS on total DCS risk as a whole. Figure 7 shows the proportion of excess DCS at ten years after HL, for male patients diagnosed in 1985 and 2000, as a function of age at diagnosis. We could alternatively quantify such a proportion over the other timescales, i.e. as a function of time since diagnosis, or indeed year at diagnosis.

Example code to conduct the previous analysis is included in the Appendix. Extended predictions such as length of stay in each state, the probability of ever visiting a state, standardised predictions and further contrasts are also provided.

Research on late effects of treatment among cancer survivors is becoming increasingly important as more patients are cured from their disease. Studying life span trajectories requires complex methods, that consider follow-up not just until the first event of interest, but through a series of intermediate states until reaching an absorbing state. Advanced multi-state methods for the cause-specific setting have been developed and are frequently used in applied research, but are unable to separate between the probability of entering an intermediate state associated with a diagnosis of cancer, and the probability of entering the same state not related to cancer. Utilising recently developed methods in multi-state modelling, where a simulation approach is taken together with parametric transition models, we have proposed one way to incorporate estimation of excess morbidity and mortality. This allows for a separation between the predicted excess illness probability and expected illness probability. We have further introduced the possibility for different transitions to have multiple time scales within a relative survival setting.

The developments provided in this paper were motivated by the study of excess morbidity and mortality of DCS among HL patients. The excess DCS mortality has been shown to decline in Sweden since the mid-1970s[14], and is no longer a major cause of death among this group of patients. Still, the risk of being diagnosed with excess DCS persists in the new millennium among HL patients in Sweden[11]. Although a reduction in excess DCS mortality does not necessarily imply a reduction in excess incidence, it remains difficult to draw conclusions on the relationship between HL and its treatment and risk of DCS together with fatality of DCS. By incorporating morbidity and mortality in the same setting, and separating between the incidence attributable to the cancer and that expected in the absence of cancer, we have illustrated one way to overcome this issue. Our approach has been implemented in the merlin and multistate packages in Stata, and code illustrating how to produce the results presented in this paper, can be found in the Appendix.

For illustrative purposes, we chose to simplify the models in our example. Only main effects models were considered, and no information on clinical covariates were included. Therefore, caution should be taken in terms of interpreting the results presented herein. It is important to point out that more complex, and realistic, models could be fitted to each of the transitions. This includes, but is not limited to, interactions between covariates and relaxing the proportional hazards assumptions [28]. As with any standard parametric modelling process, diagnostics and goodness-of-fit inference based on Wald tests, likelihood ratio tests, information criteria (AIC, BIC) can be applied and used to guide model selection for each of the transitions. Lee et al. (2017)[21] recently highlighted that standard Martingales must be adapted for the cause-specific timescale setting, and as such, further research is required to adapt them to our more general multiple timescale setting. Indeed, the example model specified for the excess incidence rate (Equation 3) here only depended on baseline age and calendar year, whereas multiple timescales could be incorporated also for this model. Moreover, the final predicted probabilities rely on four different models so in all applications it is important to assess the goodness of fit of each model. We exemplify this in the Appendix by varying the degrees of freedom; 3-5 degrees of freedom for transition 2 and 3, and 3-4 degrees of freedom for transition 4 (Supplementary Figures  8 and  9). In our illustrative example we partitioned the illness state into an expected and excess component and focused on one of several possible late effects. However, partitioning the absorbing states can be done in a similar manner given the appropriate expected mortality table, and additional illness states, e.g., secondary malignancies can be added to the model.

In terms of quantifying model results, we have given examples of how to predict transition probabilities, difference in probabilities, and proportion of excess illness over total probability of illness. Our approach can be used to predict other quantities, such as length of stay in each state (including differences and ratios thereof between patients with different covariate patterns). All of the described extensions can be implemented with the associated software package.

Several assumptions are made when using a relative survival framework. Patients are assumed to be exchangeable with the general population, conditioned on the stratification variables. If this assumption is not fulfilled, population incidence or mortality tables require additional stratification. If this is not possible, information on a control population can be used to adjust the expected rate [29]. Additionally, if we wish to interpret the excess rate as treatment-related, then we must assume that the cancer is not directly associated with illness.

Cancer survivors are at risk of a range of diseases for the remainder of their lives. To support this clinical challenge with scientific real-world evidence, we believe that our approach can be useful for predicting survivors trajectories on a population-based level and to illustrate which of different personal and clinical characteristics are at highest risk of experiencing at different time points following their cancer diagnosis.