A multistate additive relative survival semi-Markov model

Abstract

Medical researchers are often interested to investigate the relationship between explicative variables and times-to-events such as disease progression or death. Such multiple times-to-events can be studied using multistate models. For chronic diseases, it may be relevant to consider semi-Markov multistate models because the transition intensities between two clinical states more likely depend on the time already spent in the current state than on the chronological time. When the cause of death for a patient is unavailable or not totally attributable to the disease, it is not possible to specifically study the associations with the excess mortality related to the disease. Relative survival analysis allows an estimate of the net survival in the hypothetical situation where the disease would be the only possible cause of death. In this paper, we propose a semi-Markov additive relative survival (SMRS) model that combines the multistate and the relative survival approaches. The usefulness of the SMRS model is illustrated by two applications with data from a French cohort of kidney transplant recipients. Using simulated data, we also highlight the effectiveness of the SMRS model: the results tend to those obtained if the different causes of death are known.

Keywords

Disease progression kidney transplantation multistate models relative survival semi-Markov process survival analysis

1 Introduction

The method most widely used to investigate the relationship between explicative variables and times-to-events is the Cox proportional hazards model.¹ However, only a single time-to-event can be studied, which may not be relevant in a medical context. As an example, if the event under study is the progression to an advanced stage of a chronic disease, it is often difficult to consider death as non-informative censoring. Multistate models allows multiple times to event data to be studied and have been well-established for years.^2–4

Another advantage of multistate approaches is the possibility to separately study the death according to their causes. To personalize therapeutic strategies, it is important to study more specifically the treatment effects on time-to-death related to the disease under consideration. However, deaths are often studied indifferently to their cause, the cause of death being usually unavailable or hardly attributable to a unique cause. For instance, in renal transplantation, a death related to a post-transplantation cancer may be attributed to the renal chronic disease since immunosuppressive drugs represent risk factors for cancer⁵ and/or may be due to other factors independent from transplantation, such as tobacco use or hereditary factors. Relative survival (RS) models with additive hazards have been introduced as a methodological answer to this issue:^6–8 the excess mortality related to a disease is estimated by subtracting the expected mortality rates from general population life tables to the observed mortality rates. This decomposition allows an estimate of the net survival, i.e. the survival in the situation where the disease under study would be the only possible cause of death.

In light of this, it is interesting to combine both approaches, i.e. multistate model with the concept of net survival. Exploring this idea, Belot et al.⁹ considered the excess mortality related to colon cancer with competing risks data. More recently, in a general multistate framework, Huszti et al.¹⁰ proposed a time-inhomogeneous Markov relative survival (MRS) model with piecewise constant intensities. The authors demonstrated the usefulness of this model, especially when the effects of risk factors differ substantially between disease-specific mortality and expected mortality in the general population. Time-inhomogeneous Markov models consider that the transition intensities depend on the chronological time, i.e. the time since the baseline of the study. However, this assumption may not be relevant in many medical applications, for instance in chronic diseases for which time-to-death more likely depends on the time since the last disease progression than on the time since the initial disease diagnosis. It has already been demonstrated that the semi-Markovian assumption, that considers the transition intensities between two states depend on the time spent in the current state,^11,12 was particularly relevant in renal transplantation^13,14 or in AIDS.^15,16

In that context, we propose a relative survival model to describe the multistate progression of a disease considering the time spent in the current state. We call this approach SMRS (semi-Markov additive relative survival). Section 2 presents the SMRS model. Section 3 presents two applications related to the long-term evolution of kidney transplant recipients (KTRs) and describes differences with the MRS model proposed by Huszti et al.¹⁰ Finally, Section 4 offers discussions and conclusions. Of note, simulations for validation of the SMRS model properties are detailed in the online supplementary materials.

2 The SMRS model

2.1 Notation

We adopt the notation T for the chronological time from baseline and S for the duration (or sojourn time) in a state. Let $X$ the finite space of the possible clinical states. The stochastic process under consideration is ${X_{m}, T_{m}, m \in ℕ}$ , where X_m is the state of the patient after the mth transition occurring at time T_m with $T_{0} < T_{1} < \dots < T_{m}$ (T₀ = 0 and X₀ = 1 by convention). Let ε be the set of possible transitions ij with $(i, j) \in (X, X)$ , where i represents a transient state with j distinct from i. We denote by Z the overall vector of patient characteristics and Z_ij the subvector of characteristics specifically associated with the transition ij. A possible multistate representation is presented in Figure 1 (solid lines). There are two transient states (X = 1 and X = 2) and two absorbing states (X = 3 and X = O, the latter for observed death regardless of cause). We assume that no backward transition is clinically possible.

Figure 1.

The SMRS representation for a multistate model including the death related to the disease (X = E) and the death not related to the disease (X = P). Arrows for the transitions to X = E and X = P are represented with dashed lines since the two states cannot be distinguished individually.

2.2 The semi-Markov model

The semi-Markov (SM) model considers that the transition intensities between two states depend on the time already spent in the current state. In this paper, we focus on a time-homogeneous SM model. Unlike the time-inhomogeneous SM model, the time-homogeneous SM model considers that the transition intensities between two states do not depend on the time since baseline. The instantaneous hazard function specific from state X_m = i to the state X_m_+ 1 = j after a duration s, given patient characteristics Z_ij = z_ij is therefore defined by

λ_{ij} (s | z_{ij}) = {lim}_{Δ s \to 0 +} P (s \leq T_{m + 1} - T_{m} < s + Δ s, X_{m + 1} = j | T_{m + 1} - T_{m} > s, X_{m} = i, z_{ij}) / Δ s

(1)

with

Λ_{ij} (s | z_{ij}) = \int_{0}^{s} λ_{ij} (u | z_{ij}) 1 mu d u

the corresponding cumulative hazard function.

2.3 The proposed SMRS model

In additive relative survival analysis, the observed mortality hazard λ_iO(.) is assumed to be the sum of two components: the expected population mortality and the excess mortality related to the disease under study.⁸ Let X = E the death related to the disease and X = P the death related to other causes. Figure 1 shows the SMRS representation for a multistate model: whereas the usual SM model considers four observed states (solid lines), the SMRS model considers five states (dashed lines) with two types of death, although these two states cannot be distinguished individually.

The expected population mortality hazard λ_P(.) can be obtained from the death rates provided by life tables of statistical national institutes. These tables indicate the proportion of people dead in a calendar year stratified by birthdate and gender. Other strata can be included, such as ethnicity, but such additional strata are rarely available. Given the completeness of deaths from exhaustive life tables, λ_P(.) is always considered as a fixed parameter. This assumption has to be relaxed when λ_P(.) is obtained from a survival model estimated from a non-exhaustive database. In this case, sample-to-sample fluctuations should be considered.¹⁷ Let A the random variable for patient age at death, a_i his/her age observed at entry in state i, y his/her birthyear, and g his/her gender. The instantaneous hazard function for the expected mortality (not related to the disease) after a duration s in the state i is

λ_{P} (s + a_{i} | y, g) = {lim}_{Δ s \to 0 +} P (s + a_{i} \leq A < s + a_{i} + Δ s, X = P | A > s + a_{i}, y, g) / Δ s

(2)

The instantaneous hazard function of the observed mortality can be rewritten

λ_{iO} (s | z_{iE}, a_{i}, y, g) = λ_{iE} (s | z_{iE}) + λ_{P} (s + a_{i} | y, g)

(3)

The corresponding cumulative hazard function is therefore

Λ_{iO} (s | z_{iE}, a_{i}, y, g) = Λ_{iE} (s | z_{iE}) + Λ_{P} (s + a_{i} | y, g) - Λ_{P} (a_{i} | y, g)

(4)

where

Λ_{P} (s + a_{i} | y, g) - Λ_{P} (a_{i} | y, g)

represents the cumulative hazard of death between age a_i and a_i + s in the general population. The probability for a patient to sejourn at least a duration s in state i can be defined as

S_{i .} (s | z, a_{i}, y, g) = exp [- (\sum_{\binom{j : ij \in ε}{j \neq death}} Λ_{ij} (s | z_{ij})) - Λ_{iE} (s | z_{iE}) - Λ_{P} (s + a_{i} | y, g) + Λ_{P} (a_{i} | y, g)]

(5)

Based on equations (1), (3), and (5), the density function specific to transition ij, after a duration s is

f_{ij} (s | z, a_{i}, y, g) = (11_{{j \neq death}} λ_{ij} (s | z_{ij}) + 1_{{j = death}} λ_{iO} (s | z_{iE}, a_{i}, y, g)) S_{i .} (s | z, a_{i}, y, g)

(6)

The contribution of a patient to the likelihood of the SMRS model is written differently according to the last observed state. Let s_ij the duration time in state i before transition to state j, and δ_ij = 1 if the transition ij is observed, and δ_ij = 0 otherwise. For a patient in an absorbing state at his/her last time of follow-up, the contribution is

\underset{ij \in ε}{Π} {f_{ij} (s_{ij} | z, a_{i}, y, g)} δ_{ij}

. For a patient censored in the transient state k (for a duration s_k) at his/her last time of follow-up, the contribution is

\underset{ij \in ε}{Π} {f_{ij} (s_{ij} | z, a_{i}, y, g)} δ_{ij} S_{k .} (s_{k} | z, a_{k}, y, g)

The instantaneous hazard functions (including the excess death hazard λ_iE) can be modeled classically using parametric or semi-parametric distributions, and proportional or non-proportional hazard models.¹⁸ In this paper, the methods discussed are limited to parametric proportional hazard models with time-fixed covariates. The hazard ratios for transition ij (HR _ij ), calculated as the exponential of the regression coefficients, are cause-specific to the transition ij: the individuals who suffer an event k ≠ j are removed from the risk set to do a competing transition after this time event (right conditioning in equation (1)). The HR _ij is a measure of association that does not necessarily directly translate into a measure of risk without the assumption of independence between the competing events. It can be interpreted as the relative change in the cause-specific hazard for the transition ij corresponding to a 1-unit increase in the corresponding covariate.^3,19

We used the R statistical software²⁰ version 3.0.1 with optim() function to maximize the likelihood function and to compute the corresponding Hessian matrix (Nelder and Mead algorithms²¹). The complete methodology has been implemented in an R package Multistate available at http://www.divat.fr and is additionally available from the authors upon request. It includes functions for estimating SM, SMRS, Markov, and MRS models.

2.4 Comparison between observed and expected numbers of transitions

To assess the fit of the SRMS model, we proposed an ad hoc comparison between observed and expected numbers of transitions ij according to chronological time intervals $[t_{l - 1}, t_{l}] $ (l = 1, \dots, L)$ . To calculate the expected numbers of transitions between t_l₋₁ and t_l, we needed to determine for each subject the chronological time during which he/she was at risk of events. We simulated the time of dropout v using a Weibull distribution, estimated from the entire sample (the time-to-dropout being the random variable under interest and considering the time-to-absorbing event as right-censoring). Let v_h this simulated time for the subject h. The expected numbers of transitions 1j between t_l₋₁ and t_l are

e_{l, 1 j} = \sum_{h} {11_{{v_{h} > t_{l - 1}}} \int_{t_{l - 1}}^{min (v_{h}, t_{l})} f_{1 j} (u | z h, a_{1}^{h}, y h, g h) 1 mu d u}

(7)

The expected numbers of transitions jk (j ≠ 1) between t_l₋₁ and t_l are

e_{l, 1 jk} = \sum_{h} {11_{{v_{h} > t_{l - 1}}} \int_{t_{l - 1}}^{min (v_{h}, t_{l})} (\int_{0}^{v} f_{1 j} (u | z h, a_{1}^{h}, y h, g h) f_{jk} (v - u | z h, a_{j}^{h}, y h, g h) 1 mu d u) 1 mu d v}

(8)

Details of the calculations are given in Section A of the online supplementary materials.

In parallel, in order to evaluate the suitability of the parametric assumptions, we also compared for each transition the estimations of the cumulative incidence functions (CIFs) from the SMRS model and from the non-parametric Aalen–Johansen estimator.²² Estimations of the CIFs from the SMRS model for some specific strata were calculated using the mean of individuals’ CIFs, i.e. the integrals of the density functions (6). Effect of covariate on CIF should be interpreted with caution because a decrease in the CIF could be either due to a physiological effect of the exposure or to an increase in the probability of the competing failure causes.²³

3 Applications in kidney transplantation

Kidney transplantation is the treatment of choice for most patients with end-stage renal disease because survival and quality of life are often better than long-term dialysis. However, with the shortage of donor organs and the risk associated with the surgical procedure, it is of particular interest to study the risk of death related to the kidney transplant status. In particular, it is important to evaluate whether transplantation in elderly patients is associated with a dramatic increase in the mortality compared with the general population. Foucher et al.⁷ proposed, for the first time in a transplantation setting, the use of an additive regression model for relative survival. They were not able to demonstrate a significant excess of mortality related to ageing compared with the general population. They also suggested that preventing delayed graft function (DGF) may be more effective for decreasing the risk of death specifically related to the disease. Nevertheless, they only considered the time-to-death by assuming the non-informative right-censoring of returns in dialysis. Moreover, they ignored intermediate stages of disease evolution, such as the acute rejection episode (ARE) or cancer. The ARE is a common complication in the first year post-transplantation partially due to DGF. Cancer is the third cause of death in KTRs, with mortality rates higher in youngest KTR compared to the general population because of the long-term use of immunosuppressive drugs.²⁴

To illustrate the usefulness of the SMRS model, we propose studying the relationship between the DGF and the multistate evolution of KTR (Figure 1): the states X = 3 and X = O are respectively the definitive return to dialysis and the death with a functional graft, and the transient state (X = 2) is the ARE in the first application and cancer in the second application. Both applications were performed using data from DIVAT (Données Informatisées et VAlidées en Transplantation, http://www.divat.fr), a French multicentric prospective cohort of patients from six centers, older than 18 years, who have received renal transplantation and have given their informed consent for their data collection. We included 5943 adults receiving a kidney transplant between January 2000 and September 2012, without simultaneous transplanted organ. Mean recipient age was 50.9 years old (standard deviation [SD] 13.5), 3622 (61%) patients were male recipients and 1817 (31%) had a DGF. In October 2012, 625 (10.5%) patients were returned to dialysis and 351 (5.9%) were dead with a functional graft.

For each application, we compared the proposed SMRS model with two other models: the MRS model with piecewise constant intensities proposed by Huszti et al.¹⁰ and the traditional SM model without distinction between deaths (four-state SM model). The expected mortality rates of the general French population, available by age, gender and calendar year in the R package ROCt (http://www.divat.fr), were used in relative survival models. Original data and updates can also be downloaded from the Human Mortality Database (http://www.mortality.org). Three covariates were considered for each transition intensities: recipient age, recipient gender, and the DGF. The number of intervals and cutpoints related to the piecewise baseline hazard functions were determined graphically for each transition using the CIFs from the non-parametric Aalen–Johansen estimator as a reference.

3.1 The ARE as a transient state

There were 3769 (63%) living patients with a functional graft without ARE at the end of their follow-up. A total of 1198 (20%) patients had an ARE without observing death or returning to dialysis, 285 (5%) patients returned to dialysis after ARE, and 80 (1%) patients died with a functional graft after ARE. Table 1 presents the results obtained from the SMRS, MRS, and SM models. Among the three models, the SMRS model had the smallest Akaike information criterion (AIC).

Table 1.

Comparison of the multivariable SMRS, MRS and four-state SM models for the analysis of the relationship between clinical factors and the time-to-events on renal transplant data with ARE as a transient state (N = 5943). For the HRs related to age, the reference group corresponds to patients older than 65 years.

		SMRS model¹			MRS model²			Four-state SM model³
Transition	Coefficient	HR	[95% CI]	p-value	HR	[95% CI]	p-value	HR	[95% CI]	p-value
	Age < 35 years	1.78	[1.44–2.21]	0.0000	1.79	[1.44–2.22]	0.0000	1.78	[1.43–2.20]	0.0000
12	Age 35–50 years	1.55	[1.28–1.89]	0.0000	1.56	[1.28–1.90]	0.0000	1.55	[1.27–1.88]	0.0000
(transplantation	Age 50–65 years	1.28	[1.05–1.55]	0.0139	1.28	[1.05–1.56]	0.0131	1.27	[1.05–1.55]	0.0153
to ARE ⁴ )	Male recipient	1.15	[1.04–1.27]	0.0072	1.15	[1.04–1.27]	0.0080	1.15	[1.05–1.27]	0.0071
	Delayed graft function	1.23	[1.12–1.36]	0.0001	1.23	[1.12–1.36]	0.0001	1.23	[1.12–1.36]	0.0001
	Age < 35 years	0.45	[0.31–0.65]	0.0000	0.45	[0.31–0.65]	0.0000	0.44	[0.31–0.64]	0.0000
13	Age 35–50 years	0.47	[0.34–0.65]	0.0000	0.47	[0.34–0.65]	0.0000	0.47	[0.33–0.65]	0.0000
(transplantation to	Age 50–65 years	0.53	[0.40–0.72]	0.0000	0.53	[0.40–0.72]	0.0000	0.53	[0.40–0.71]	0.0000
return to dialysis)	Male recipient	0.87	[0.70–1.07]	0.1894	0.86	[0.70–1.07]	0.1814	0.86	[0.70–1.07]	0.1842
	Delayed graft function	2.47	[1.99–3.06]	0.0000	2.46	[1.98–3.05]	0.0000	2.46	[1.98–3.05]	0.0000
	Age < 35 years	0.06	[0.01–0.31]	0.0009	0.06	[0.01–0.30]	0.0007	0.06	[0.03–0.13]	0.0000
1E/1O	Age 35–50 years	0.34	[0.19–0.61]	0.0003	0.34	[0.19–0.61]	0.0003	0.22	[0.15–0.31]	0.0000
(transplantation	Age 50–65 years	0.55	[0.32–0.93]	0.0224	0.55	[0.32–0.93]	0.0233	0.45	[0.33–0.60]	0.0000
to death ⁵ )	Male recipient	0.78	[0.53–1.16]	0.2304	0.78	[0.53–1.16]	0.2297	1.14	[0.89–1.48]	0.2894
	Delayed graft function	3.02	[1.96–4.64]	0.0000	3.02	[1.96–4.65]	0.0000	1.93	[1.52–2.44]	0.0000
	Age < 35 years	0.88	[0.56–1.39]	0.5851	0.87	[0.55–1.36]	0.5402	0.88	[0.56–1.38]	0.5727
23	Age 35–50 years	0.71	[0.46–1.10]	0.1284	0.71	[0.46–1.09]	0.1223	0.71	[0.46–1.09]	0.1234
(ARE ⁴ to return	Age 50–65 years	0.64	[0.42–0.99]	0.0447	0.64	[0.42–0.99]	0.0480	0.64	[0.41–0.98]	0.0424
to dialysis)	Male recipient	0.64	[0.51–0.82]	0.0002	0.64	[0.50–0.80]	0.0002	0.64	[0.51–0.81]	0.0002
	Delayed graft function	1.61	[1.28–2.04]	0.0001	1.62	[1.28–2.04]	0.0001	1.61	[1.27–2.04]	0.0001
	Age < 35 years	0.18	[0.04–0.79]	0.0221	0.17	[0.04–0.69]	0.0128	0.11	[0.04–0.30]	0.0000
2E/2O	Age 35–50 years	0.36	[0.11–1.18]	0.0913	0.31	[0.10–0.98]	0.0480	0.29	[0.14–0.58]	0.0005
(ARE ⁴ to death ⁶ )	Age 50–65 years	0.67	[0.22–2.02]	0.4757	0.61	[0.22–1.70]	0.3529	0.55	[0.29–1.02]	0.0569
	Male recipient	0.44	[0.22–0.89]	0.0216	0.47	[0.24–0.94]	0.0294	0.70	[0.44–1.10]	0.1128
	Delayed graft function	2.49	[1.23–5.04]	0.0119	2.57	[1.29–5.10]	0.0077	1.88	[1.20–2.95]	0.0052
Log-likelihood		−7715.8			−7734.3			−7727.2
AIC		15,521.6			15,562.5			15,544.3

Cutpoints for baseline hazards $(s_{l}, l = 1, \dots, L)$ and log(constant) $log (constant) (log (a_{l}), l = 1, \dots, L)$ are: transition 12 (0, 20 days, 100 days, 375 days, 4 years, 8 years), (−0.07,−1.15,−2.70,−4.28,−4.82,−5.91); transition 13 (0, 20 days, 120 days, 3 years), (−1.60,−3.14,−4.15,−3.55); transition 1E (0, 20 days, 6 months, 5 years), (−2.90,−3.81,−4.59,−4.05); transition 23 (0, 20 days, 1 year), (−0.83,−2.49,−2.77); transition 2E (0, 90 days, 3 years), (−3.67,−4.57,−3.52). ² Cutpoints for baseline hazards $(t_{m}, m = 1, \dots, M)$ and log(constant) $log (constant) (log (b_{m}), m = 1, \dots, M)$ are: transition 12 (0, 20 days, 100 days, 375 days, 4 years, 8 years), (−0.07,−1.15,−2.70,−4.28,−4.82,−5.89); transition 13 (0, 20 days, 120 days, 3 years), (−1.60,−3.13,−4.15,−3.54), transition 1E (0, 20 days, 6 months, 5 years), (−2.89,−3.81,−4.60,−4.05); transition 23 (0, 20 days), (−0.59,−2.62); transition 2E (0, 20 days, 100 days, 4 years, 8 years), (−10.71,−3.21,−6.00,−4.27,−3.21,−4.41). ³ Cutpoints for baseline hazards $(s_{l}, l = 1, \dots, L)$ and log(constant) $(log (a_{l}), l = 1, \dots, L)$ are: transition 12 (0, 20 days, 100 days, 375 days, 4 years, 8 years), (−0.07,−1.14,−2.70,−4.27,−4.82,−5.91); transition 13 (0, 20 days, 120 days, 3 years), (−1.60,−3.14,−4.14,−3.54); transition 1O (0, 20 days, 6 months, 5 years), (−2.58,−3.29,−3.67,−3.29); transition 23 (0, 20 days, 1 year), (−0.83,−2.49,−2.77); transition 2O (0, 90 days, 3 years), (−3.12,−3.79,−3.02). ⁴ ARE: acute rejection episode ⁵ Transplantation to death related to the disease (1E) for the SMRS and MRS models, transplantation to observed death (1O) for the four−state SM model. ⁶ ARE to death related to the disease (2E) for the SMRS and MRS models, ARE to observed death (2O) for the four−state model.

Results from the SMRS and MRS models were similar for transitions from the first state. This was expected as duration and chronological times coincide for this transition. Nevertheless the results for these two models were also close from the transitions to ARE (X = 2). This may be also explained by the short time to ARE, leading to close duration and chronological times after ARE.

As expected, estimations for transitions into clinical events other than death were quite similar between the SMRS model and the four-state SM model. More interestingly, the estimations from the SMRS model differed from the SM model for transitions to death. This was expected because some deaths are related to the kidney transplant status. The effect of ageing on death related to the transplantation (SMRS model) was lower than those for all-cause death (four-state SM model) but remained important indicating that the role of ageing observed in the SM model might not only be explained by the expected effect of age as in the general population. For example, in the SMRS model, the HRs of death related to the disease without ARE, varied from 0.06 (95% confidence interval (CI) [0.01–0.31]) for recipients younger than 35 years old to 0.55 (95% CI [0.32–0.93]) for recipients between 50 and 65 years old compared with recipients older than 65 years of age. In the SM model these HR varied from 0.06 (95% CI [0.03–0.13]) to 0.45 (95% CI [0.33–0.60]). Similar trends could be underlined for the transition from ARE to death. In contrast, the associations between DGF and the different risks of transition into death were greater in the SMRS than in the SM model. In the SMRS model, the corresponding HR for death related to the disease was estimated to be 3.02 (95% CI [1.96–4.64]) without ARE, and to be 2.49 (95% CI [1.23–5.04]) after ARE. In the SM model, these HRs were respectively estimated to be 1.93 (95% CI [1.52–2.44]) and 1.88 (95% CI [1.20–2.95]). The differences between both approaches are consistent since the DGF cannot be a risk factor of death related to a cause other than the disease. Finally, we observed that women seemed to have a higher risk than men to die due to the disease, especially after ARE, when the association was embedded in the all-cause model (the mortality for men being higher than for women in general population).

The CIFs derived from the SMRS model and from the non-parametric Aalen–Johansen estimator are presented in Section B of the online supplementary materials. Regardless of the type of event or group of patients, the SMRS model estimates tended to fit the non-parametric ones. Therefore, the parametric assumptions for durations and the proportionality of hazards seemed reasonable.

To assess the homogeneity assumption of the SMRS model regarding chronological time, we compared the observed and expected numbers of transitions using five intervals of chronological times defined with the quantiles of the observed times of transitions (Table 2). The differences for transitions from state X = 1 were fully acceptable with either relative differences lower than 10% or absolute differences lower than 10 transitions. The SMRS model tended to predict ARE events too early: the slightly excessive expected number of ARE before 3 months ( + 8.7%) was offset by a deficit after this period. The expected numbers of returns to dialysis and deaths, without ARE were close to clinical observations. Yet the relative differences for transitions after ARE were large for some time intervals but there were few transitions. The larger errors observed for returns to dialysis after ARE can also be a consequence of the earlier ARE predicted by the model. There was no evidence of increased differences between observed and expected numbers of transitions with the chronological time, indicating that the homogeneity assumption of the SMRS model regarding the chronological time seemed reasonable.

Table 2.

Contingency table for the observed and expected numbers of transitions derived from the SMRS model with ARE as a transient state (N = 5943).

Chronological		Transition
time (months)	Counts¹	12	13	14	23	24
[0.0;0.6]	Observed	492	46	13	4	0
	Expected	486.3	46.7	13.1	4.4	0.3
	Absolute difference	−5.7	0.7	0.1	0.4	0.3
	Relative difference	−1.2%	1.4%	1.1%	10.0%	ND²
[0.7;3.3]	Observed	477	30	23	7	4
	Expected	518.4	31.5	19.8	17.0	2.5
	Absolute difference	41.4	1.5	−3.2	10.0	−1.5
	Relative difference	8.7%	5.1%	−13.4%	143.2%	−37.3%
[3.3;12.5]	Observed	361	47	43	45	7
	Expected	352.3	40.7	45.5	50.5	9.1
	Absolute difference	−8.7	−6.3	2.5	5.5	2.1
	Relative difference	−2.4%	−13.4%	5.9%	12.3%	29.8%
[12.6;45.5]	Observed	163	86	88	118	28
	Expected	149.6	88.3	84.7	100.5	26.5
	Absolute difference	−13.4	2.3	−3.3	−17.5	−1.5
	Relative difference	−8.2%	2.7%	−3.7%	−14.8%	−5.3%
[45.6;148.7]	Observed	70	131	104	111	41
	Expected	68.6	119.3	106.1	98.1	46.4
	Absolute difference	−1.4	−11.7	2.1	−12.9	5.4
	Relative difference	−2.0%	−8.9%	2.0%	−11.6%	13.1%

¹Absolute difference = Expected counts – Observed counts

Relative difference = (Expected counts – Observed counts)/Observed counts; ² Not defined.

3.2 The post-transplant cancer as a transient state

There were 4537 (76%) living patients with a functional graft without cancer at the end of their follow-up. A total of 554 (9%) patients had cancer without observing death or returning to dialysis, 40 (1%) patients returned to dialysis after cancer, and 84 (1%) patients died with a functional graft after cancer. Table 3 presents the results obtained from the SMRS, MRS, and SM models. Among the three models, the SMRS model had the smallest AIC.

Table 3.

Comparison of the multivariable SMRS, four-state SM and MRS models for the analysis of the relationship between clinical factors and the time-to-events on renal transplant data with cancer as a transient state (N = 5943). For the HRs related to age, the reference group corresponds to patients lower than 55 years.

		SMRS model¹			MRS model²
Transition	Coefficient	HR	[95% CI]	p-value	HR	[95% CI]	p-value	HR	[95% CI]	p-value
12	Age > 55 years	3.27	[2.74–3.90]	0.0000	3.23	[2.70–3.85]	0.0000	3.26	[2.74–3.89]	0.0000
(transplantation	Male recipient	1.38	[1.15–1.64]	0.0005	1.35	[1.13–1.61]	0.0009	1.38	[1.15–1.64]	0.0005
to cancer)	Delayed Graft Function	1.10	[0.92–1.31]	0.3006	1.10	[0.92–1.31]	0.2986	1.10	[0.92–1.31]	0.2987
13	Age > 55 years	1.24	[1.04–1.48]	0.0129	1.26	[1.05–1.50]	0.0081	1.24	[1.04–1.48]	0.0131
(transplantation to	Male recipient	0.80	[0.69–0.94]	0.0086	0.81	[0.69–0.94]	0.0101	0.80	[0.69–0.94]	0.0084
return to dialysis)	Delayed Graft Function	2.10	[1.80–2.46]	0.0000	2.11	[1.80–2.47]	0.0000	2.10	[1.79–2.46]	0.0000
1E/1O	Age > 55 years	3.72	[2.37–5.84]	0.0000	3.48	[2.31–5.25]	0.0000	3.65	[2.83–4.71]	0.0000
(transplantation	Male recipient	0.55	[0.36–0.84]	0.0052	0.56	[0.38–0.83]	0.0042	0.91	[0.72–1.15]	0.4578
to death)	Delayed graft function	3.55	[2.26–5.58]	0.0000	3.35	[2.22–5.05]	0.0000	2.28	[1.80–2.88]	0.0000
23	Age > 55 years	1.17	[0.61–2.23]	0.6382	1.42	[0.75–2.72]	0.2808	1.17	[0.61–2.23]	0.6386
(cancer to	Male recipient	0.97	[0.50–1.90]	0.9399	1.20	[0.62–2.34]	0.5927	0.97	[0.50–1.89]	0.9345
return to dialysis)	Delayed graft function	2.75	[1.47–5.15]	0.0016	3.14	[1.71–5.77]	0.0002	2.76	[1.47–5.17]	0.0016
2E/2O	Age > 55 years	1.05	[0.61–1.82]	0.8555	1.14	[0.67–1.94]	0.6260	1.32	[0.84–2.07]	0.2360
(cancer to death)	Male recipient	1.18	[0.65–2.12]	0.5868	1.50	[0.82–2.76]	0.1841	1.32	[0.81–2.15]	0.2723
	Delayed graft function	1.61	[0.95–2.73]	0.0821	1.86	[1.12–3.09]	0.0168	1.53	[0.99–2.35]	0.0537
Log-likelihood		−6907.4			−6920.2			−6922.5
AIC		13,892.8			13,910.5			13,923.0

¹Cutpoints for baseline hazards $(s_{l}, l = 1, \dots, L)$ and log(constant) ( $log (a_{l}), l = 1, \dots, L$ ) are: transition 12 (0, 6 months, 1 year, 4 years, 8 years), (−4.67,−4.67,−4.53,−4.23,−4.07); transition 13 (0, 20 days, 120 days, 3 years), (−2.10,−3.69,−4.13,−3.79); transition 1O (0, 20 days, 6 months, 5 years), (−4.37,−5.34,−6.30,−5.70); transition 23 (0, 30 days, 60 days, 6 months, 1 year, 5 years), (−3.63,−2.66,−4.50,−4.51,−4.30,−3.68); transition 2O (0, 90 days, 1 year, 2 years, 4 years), (−2.07,−3.30,−3.81,−4.01,−3.93). ²Cutpoints for baseline hazards ( $t_{m}, m = 1, \dots, M$ ) and log(constant) ( $log (b_{m}), m = 1, \dots, M$ ) are: transition 12 (0, 6 months, 1 year, 4 years, 8 years), (−4.63,−4.61,−4.52,−4.21,−3.97); transition 13 (0, 20 days, 120 days, 3 years), (−2.08,−3.67,−4.15,−3.81); transition 1E (0, 20 days, 6 months, 5 years), (−4.12,−5.18,−6.19,−5.53); transition 23 (0, 6 months, 2 years, 3 years), (−2.76,−3.92,−5.82,−4.42); transition 2E (0, 1 year, 8 years), (−2.31,−3.64,−4.47). ³Cutpoints for baseline hazards ( $s_{l}, l = 1, \dots, L$ ) and log(constant) ( $log (a_{l}), l = 1, \dots, L$ ) are: transition 12 (0, 6 months, 1 year, 4 years, 8 years), (−4.67,−4.67,−4.53,−4.22,−4.07); transition 13 (0, 20 days, 120 days, 3 years), (−2.10,−3.69,−4.12,−3.79); transition 1O (0, 20 days, 6 months, 5 years), (−4.23,−4.98,−5.52,−5.06), transition 23 (0, 30 days, 60 days, 6 months, 1 year, 5 years), (−3.63,−2.66,−4.50,−4.52,−4.30,−3.68); transition 2O (0, 90 days, 1 year, 2 years, 4 years), (−2.22,−3.29,−3.79,−3.77,−3.68). ⁴ Transplantation to death related to the disease (1E) for the SMRS and MRS models, transplantation to observed death (1O) for the four-state SM model. ⁵ Cancer to death related to the disease (2E) for the SMRS and MRS models, cancer to observed death (2O) for the four-state model.

As in the first application, results obtained with the SMRS and MRS models were close for transitions from the initial state. In contrast, regression coefficient for transitions after cancer differed importantly. For example, in the SMRS model, the estimated HRs related to DGF were 2.75 (95% CI [1.47–5.15]) for the return to dialysis, and 1.61 (95% CI [0.95–2.73]) for the death related to the disease. In the MRS model, these HRs were respectively 3.14 (95% CI [1.71–5.77]) and 1.86 (95% CI [1.12–3.09]). Additional reason to favor the SMRS model was that transition times after cancer were not significantly associated with the duration before cancer (p = 0.69 for return to dialysis and p = 0.59 for death). Similarly to the first application, differences were observed between the SMRS model and the four-state SM model for the transitions to death: (i) the regression coefficients related to DGF were greater in the SMRS than in the SM model; and (ii) when no cancer was diagnosed post-transplantation, gender appeared as an important determinant of death related to the disease whereas this factor was not significantly associated with the all-cause mortality.

The CIFs derived from the SMRS model and from the non-parametric Aalen–Johansen estimator are presented in Section C of the online supplementary materials. Similarly to the first application, the parametric assumptions for durations and the proportionality of hazards seemed reasonable. As described in Table 4, the differences between observed and expected numbers of transitions for the SMRS model were fully acceptable before 24 months with either relative differences lower than 10% or absolute differences lower than 10 transitions. The SMRS model tended to underestimate the number of events between 24 and 69 months post-transplantation and overestimate their incidence after this time, although this was reasonable.

Table 4.

Contingency table for the observed and expected numbers of transitions derived from the SMRS model with cancer as a transient state (N = 5943).

Chronological		Transition
time (months)	Counts¹	12	13	14	23	24
[0.0;6.6]	Observed	67	134	60	3	0
	Expected	65.8	122.4	59.2	0.9	2.6
	Absolute difference	−1.2	−11.6	−0.8	−2.1	2.6
	Relative difference	−1.8%	−8.7%	−1.3%	−70.0%	ND²
[6.6;24.0]	Observed	144	111	60	5	18
	Expected	140.0	120.6	61.6	4.5	12.7
	Absolute difference	−4.0	9.6	1.6	−0.5	−5.3
	Relative difference	−2.8%	8.6%	2.7%	−2.5%	−29.4%
[24.0;45.3]	Observed	125	125	48	4	18
	Expected	115.7	110.3	51.3	7.5	17.1
	Absolute difference	−9.3	−14.7	3.3	3.5	−0.9
	Relative difference	−7.4%	−11.8%	6.9%	87.5%	−5.0%
[45.3;69.0]	Observed	110	110	51	11	23
	Expected	99.1	89.8	39.5	8.8	18.5
	Absolute difference	−10.6	−20.2	−11.5	−2.2	−4.5
	Relative difference	−9.6%	−18.4%	−22.5%	−20.0%	−19.6%
[69.0;148.7]	Observed	108	105	48	17	25
	Expected	126.5	105.5	59.6	25.0	41.5
	Absolute difference	18.5	0.5	11.6	8.0	16.5
	Relative difference	17.1%	0.5%	24.2%	47.0%	66.0%

¹Absolute difference = Expected counts – Observed counts

Relative difference = (Expected counts – Observed counts)/Observed counts ² Not defined

4 Discussion

We have proposed an additive relative survival analysis using a semi-Markovian multistate model. The merging of relative survival and multistate analyses concepts may be very informative for many clinical epidemiology studies. The applications of our model to study the evolution of KTRs illustrate in practice the usefulness of the SMRS model to better describe the incidence of deaths related to the disease and the associated risk factors. Indeed, the risk/benefit of transplantation of organs into elderly patients has been debated for a long time, particularly in light of the shortage of donor organs and the higher risk associated with the surgical procedure. According to our results, recipient age remains a major risk factor of death related to the disease. Our analyses confirm the well-known consequences of the DGF after kidney transplantation²⁵ and provide additional information on the mechanistic pathways implicated. In particular, the SMRS model described a correlation between the DGF and the time-to-death related to the disease, and this was regardless of the occurrence of ARE. The fit of the SMRS model and the homogeneity assumption regarding the chronological time seemed reasonable. This supports the use of semi-Markov multistate models with KTRs data. Although the parametric assumptions for durations appeared appropriate in our application, there may be some diseases where the non-parametric distributions fit better.¹⁵

Huszti et al.¹⁰ recently developed a time-inhomogeneous MRS model. However, in some medical contexts, the transitions may be more likely associated with the time spent in the current state than to the chronological time. Our first application shows that both SMRS and MRS models tend to have similar results. This result may be explained by early ARE (X = 2), the majority occurring during the first year post-transplantation. The consequence is a close distribution between chronological and duration times. In contrast, results from the SMRS and MRS models differed in the second application, i.e. when the transient state is cancer. Cancer is a long-term complication and therefore chronological time and duration time after the occurrence of cancer differed importantly. In both applications, it turns out that the SMRS model was preferred regarding the AIC. From our point of view, the SMRS model has several benefits. First of all, the semi-Markovian assumption corresponds to a clinical reality in many chronic diseases: the risk of death more likely depends on the time since the last disease stage evolution (cancer, infection, disease recurrence, etc.) than on the time since the disease diagnosis. Second, it is straightforward to include the time spent in the first state as a covariate of transitions from X = 2. This takes the chronological time in the SMRS model into consideration. Third, the check of parametric assumptions in the SMRS model is easier than in the MRS model: the instantaneous hazard functions depend on the duration time which enables to evaluate the suitability of the parametric assumptions by comparing the CIFs estimates with the non-parametric ones.

From simulations (see the online supplementary materials), we have demonstrated that the performances of the SMRS model are as good as a SM model where the causes of death are known (five-state model), and are better than the traditional SM model (where the causes of death are not distinguished) to estimate the effects of risk factors on the mortality related to the disease. In our SMRS model, the variability of estimations associated to transitions into death was greater than for an all-cause death model. This is due to the fact that net estimations attempt to only use information regarding deaths related to the disease, whereas the all-cause estimations are based on all observed deaths. A similar increase was reported previously.^10,26 Large standard deviations could have indicated that few deaths are related to the disease under study.

It is important to note that we have used a cause-specific writing for hazard functions. The cause-specific hazard ratios calculated are a valid measure of the apparent associations between a covariate and the relative instantaneous hazard rate given that individuals have survived all competing events until time t. A marginal (or net) interpretation is only possible in the hypothetical situation that the competing event time distributions and censoring distribution are independent but this assumption cannot be verified from the observable data.^3,19,27,28

In a relative survival framework, the expected hazard of death is calculated from general population life tables. Such estimates are known to be biased since they also include mortality from the disease under study. However, the bias is sufficiently small to be insignificant for most applications and particularly for kidney transplant data since mortality among KTRs is thought to constitute only a small proportion of total mortality.²⁹ In our application, the use of general population life tables may also be discussed because the expected mortality for KTRs can also be the mortality observed in patients with renal disease who have not benefited from transplantation. To date, no lifetime table exists for such population.

As a conclusion, this study shows the usefulness of the SMRS model to investigate the relationship between covariates and a chronic disease progression, particularly toward the death related to the disease. It offers additional results to traditional multistate models, in particular the interpretation of the mechanism associated with death related to the disease. The complete methodology has been implemented in an R package to make these statistical tools more user friendly. Future developments would be to implement semi-parametric or non-parametric multistate relative survival models for cases where parametric assumptions for durations are not relevant.

Supplemental Material

sj-pdf-1-smm-10.1177_0962280215586456 - Supplemental material for A multistate additive relative survival semi-Markov model

Supplemental material, sj-pdf-1-smm-10.1177_0962280215586456 for A multistate additive relative survival semi-Markov model by Florence Gillaizeau, Etienne Dantan, Magali Giral, Yohann Foucher, Professor Marie Davidian in Statistical Methods in Medical Research

Footnotes

5 Supplementary materials

The reader is referred to the online supplementary materials (available at: ) for technical appendices, additional simulations and supplementary graphics. Section A details the calculations of expected numbers of transitions from the SMRS model (presented in Section 2.4). Sections B and C present the CIFs estimates provided by the SMRS model and the Aalen–Johansen non-parametric estimator on data from KTRs with ARE as a transient state (Section B) and cancer as a transient state (Section C) according to: DGF, recipient age, and recipient gender. Section D details the design and results of six simulation scenarios to assess the performances of the SMRS model and compare them to the five-state SM model (the causes of death are known), and the four-state SM model (all deaths are considered regardless their cause).

Acknowledgements

The authors would like to thank the Roche laboratory for partially supporting the DIVAT cohort, and participants of the DIVAT network for their contribution to ensure the maintenance of the cohort and the development of research projects.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the French Ministry of Health [grant number PHRC PROG/11/85 2011 to Pr Magali Giral].

References

Cox

Oakes

. Analysis of survival data, New York: Chapman and Hall, 1984.

Hougaard

. Multi-state models: a review. Lifetime Data Anal 1999; 5(3): 239–264.

Putter

Fiocco

Geskus

. Tutorial in biostatistics: competing risks and multi-state models. Stat Med 2007; 26(11): 2389–2430.

Meira-Machado

de Uña-Álvarez

Cadarso-Suárez

, et al. Multi-state models for the analysis of time-to-event data. Stat Methods Med Res 2009; 18(2): 195–222.

Dantal

Hourmant

Cantarovich

, et al. Effect of long-term immunosuppression in kidney-graft recipients on cancer incidence: randomised comparison of two cyclosporin regimens. Lancet 1998; 351(9103): 623–628.

Hakulinen

Tenkanen

. Regression analysis of relative survival rates. J R Stat Soc Ser C Appl Stat 1987; 36: 309–317.

Esteve

Benhamou

Croasdale

, et al. R : elements for further discussion. Stat Med 1990; 9(5): 529–538.

Perme

Stare

Estève

. On estimation in relative survival. Biometrics 2012; 68(1): 113–120.

Belot

Remontet

Launoy

, et al. Competing risk models to estimate the excess mortality and the first recurrent-event hazards. BMC Med Res Methodol 2011; 11: 78–78.

10.

Huszti

Abrahamowicz

Alioum

, et al. Relative survival multistate Markov model. Stat Med 2012; 31(3): 269–286.

11.

Lagakos

Sommer

Zelen

. Semi-Markov models for partially censored data. Biometrika 1978; 65(2): 311–317.

12.

Huzurbazar

. Multistate models, flowgraph models, and semi-Markov processes. Commun Stat Theory Methods 2004; 33(3): 457–474.

13.

Foucher

Giral

Soulillou

, et al. A semi-Markov model for multistate and interval-censored data with multiple terminal events. application in renal transplantation. Stat Med 2007; 26(30): 5381–5393.

14.

Foucher

Giral

Soulillou

, et al. A flexible semi-Markov model for interval-censored data and goodness-of-fit testing. Stat Methods Med Res 2010; 19(2): 127–145.

15.

Joly

Commenges

. A penalized likelihood approach for a progressive three-state model with censored and truncated data: application to AIDS. Biometrics 1999; 55(3): 887–890.

16.

Foucher

Mathieu

Saint-Pierre

, et al. A semi-Markov model based on generalized weibull distribution with an illustration for HIV disease. Biom J 2005; 47(6): 825–833.

17.

Trebern-Launay

Giral

Dantal

, et al. Comparison of the risk factors effects between two populations: two alternative approaches illustrated by the analysis of first and second kidney transplant recipients. BMC Med Res Methodol 2013; 13: 102–102.

18.

Giorgi

Abrahamowicz

Quantin

, et al. A relative survival regression model using B-spline functions to model non-proportional hazards. Stat Med 2003; 22(17): 2767–2784.

19.

Lau

Cole

Gange

. Competing risk regression models for epidemiologic data. Am J Epidemiol 2009; 170(2): 244–256.

20.

R Core Team . R: A language and environment for statistical computing, Vienna, Austria: R Foundation for Statistical Computing, 2013.

21.

Nelder

Mead

. A simplex-method for function minimization. Computer journal 1965; 7(4): 308–313.

22.

Aalen

Johansen

. An empirical transition matrix for non-homogeneous Markov chains based on censored observations. Scand J Statist 1978; 5(3): 141–150.

23.

Latouche

Boisson

Chevret

, et al. Misspecified regression model for the subdistribution hazard of a competing risk. Stat Med 2007; 26(5): 965–974.

24.

Hall

Engels

Pfeiffer

, et al. Association of antibody induction immunosuppression with cancer after kidney transplantation. Transplantation 2015; 99(5): 1051–1051. –7.

25.

Yarlagadda

Coca

Formica

, et al. Association between delayed graft function and allograft and patient survival: a systematic review and meta-analysis. Nephrol Dial Transplant 2009; 24(3): 1039–1047.

26.

Le Teuff

Abrahamowicz

Bolard

, et al. Comparison of Cox’s and relative survival models when estimating the effects of prognostic factors on disease-specific mortality: a simulation study under proportional excess hazards. Stat Med 2005; 24(24): 3887–3909.

27.

Tsiatis

. A nonidentifiability aspect of the problem of competing risks. Proc Natl Acad Sci USA 1975; 72(1): 20–22.

28.

Latouche

Beyersmann

Fine

. Comments on ‘Analysing and interpreting competing risk data’. Stat Med 2007; 26(19): 3676–3679.

29.

Talback

Dickman

. Estimating expected survival probabilities for relative survival analysis – exploring the impact of including cancer patient mortality in the calculations. Eur J Cancer 2011; 47(17): 2626–2632.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.58 MB