# Integrating relative survival in multi-state models – a non-parametric
approach

Damjan Manevski  Institute for Biostatistics and Medical Informatics, Faculty
of Medicine, University of Ljubljana, Vrazov trg 2, 1000 Ljubljana, Slovenia
Hein Putter  Leiden University Medical Center, Einthovenweg 20 2333 ZC Leiden,
the Netherlands  Maja Pohar Perme  Institute for Biostatistics and Medical
Informatics, Faculty of Medicine, University of Ljubljana, Vrazov trg 2, 1000
Ljubljana, Slovenia  Edouard F. Bonneville  Leiden University Medical Center,
Einthovenweg 20 2333 ZC Leiden, the Netherlands  Johannes Schetelig  Medical
Clinic I, University Hospital, Fetscherstrasse 74, 01307 Dresden, Germany
Liesbeth C. de Wreede  Leiden University Medical Center, Einthovenweg 20 2333
ZC Leiden, the Netherlands

###### Abstract

Multi-state models provide an extension of the usual survival/event-history
analysis setting. In the medical domain, multi-state models give the
possibility of further investigating intermediate events such as relapse and
remission. In this work, a further extension is proposed using relative
survival, where mortality due to population causes (i.e. non-disease-related
mortality) is evaluated. The objective is to split all mortality in disease
and non-disease-related mortality, with and without intermediate events, in
datasets where cause of death is not recorded or is uncertain. To this end,
population mortality tables are integrated into the estimation process, while
using the basic relative survival idea that the overall mortality hazard can
be written as a sum of a population and an excess part. Hence, we propose an
upgraded non-parametric approach to estimation, where population mortality is
taken into account. Precise definitions and suitable estimators are given for
both the transition hazards and probabilities. Variance estimating techniques
and confidence intervals are introduced and the behaviour of the new method is
investigated through simulations. The newly developed methodology is
illustrated by the analysis of a cohort of patients followed after an
allogeneic hematopoietic stem cell transplantation. The work is also
implemented in the R package mstate.

##  1 Introduction

Multi-state models provide a framework for simultaneously analyzing competing
events and sequences of events. In the medical field, these models help to
study the impact of intermediate events on the prognosis of patients, and
allow to estimate separate probabilities of death with and without the
intermediate event at multiple time horizons. Interest can also be in
distinguishing between deaths due to the studied disease and its treatment
(excess mortality) and deaths due to other (population) causes. This is of
particular interest in an older patient population where the risk of dying due
to other causes is high and may considerably contribute to the overall
proportion of deaths. However, cause of death is often not reported [1],
unreliable [2, 3] or cannot unequivocally be attributed to the disease or
other causes.

Our motivating example comes from the study of outcomes for patients with
myelodysplastic syndromes (MDS) or secondary acute myeloid leukemia (sAML)
after an allogeneic hematopoietic stem cell transplantation (alloHCT) [4].
Patients with these indications represent the oldest major patient population
referred to alloHCT with a median age of 58 years at alloHCT (registry data
from 2012). These patients are at significant risk for two competing failures:
relapse of the underlying disease and non-relapse mortality (NRM), which for a
large part is due to the transplantation and pre-treatment. The occurence of
relapse leads to a very poor prognosis. We investigated the contribution of
population mortality to both death after relapse and NRM. This enabled us to
estimate the probability of excess NRM, which especially for older patients
and for long-term outcomes may provide a better estimate for treatment-related
mortality than all NRM. The current paper discusses and extends the model
introduced in that paper.

When reliable information on causes of death is not available in the medical
data, one needs to address this issue using external data. This is the core
idea of the relative survival field [5]: the observed data is merged with
general population mortality tables [6] to indirectly enable the estimation of
cause-specific hazards and probabilities [7, 8]. In this work, we extend the
ideas of the relative survival methodology to a Markov time-inhomogeneous
multi-state model where we focus on non-parametric estimation of transition
hazards and probabilities. We provide a theoretical foundation, a study of
behaviour of the methods under different scenarios and a software
implementation in R that makes the methodology readily available to other
users.

These measures have also been the focus of the relative survival analysis
field, which has seen considerable development in the last decade [5, 9, 10].
With the generalisation to multi-state models, a new level of complexity is
added. The main contribution of this paper is to discuss which measures may be
of main interest in this setting and how to estimate them. We propose
estimators for cumulative hazards and transition probabilities and study their
properties, including their (asymptotic) variances. To the best of our
knowledge, this is the first extension of relative survival to non-parametric
estimation in a Markov multi-state model, whereas other semi-parametric [11,
12, 13] or parametric [14, 15] approaches have already been proposed. An
extensive simulation study of properties of the estimators in terms of bias,
standard error and coverage probability for confidence intervals has thus been
performed. The simulation study is also the first (as to our knowledge) to
extensively investigate the behaviour of different confidence interval methods
for non-parametric multi-state models.

The article is structured as follows: in Section 2 we present the main
theoretical grounds of our work. In Section 3 we present simulation results by
which this new approach is evaluated. In Section 4 we will reanalyse the
motivating dataset to illustrate the behaviour and interpretation of the new
model. This also serves as a step-by-step overview of how such an analysis can
be performed in R using the newly developed R code. In Section 5 we discuss
the main findings of the article and provide conclusions.

##  2 Extended multi-state model integrating relative survival

We start by introducing the theoretical background of multi-state models and
relative survival needed for our work in Sections 2.1 and 2.2. We then present
our proposed extension in Section 2.3.

###  2.1 General multi-state model

A multi-state model is a stochastic process $(X(t),t\in[0,\tau))$ with a
finite state space $S$ of size $K$, $\tau\in{\rm I\\!R}^{+}$, where the sample
paths are assumed to be right-continuous, i.e. $X(t+)=X(t)$. Thus, $X(t)$
represents the state at which the process (or a given individual) is present
at time $t$. A transition from state $h$ to state $j$ takes place when an
individual at risk in state $h$ experiences an event that makes them enter
state $j$.

By assuming $X(t)$ is a Markov process we can define transition probabilities
in the following way

$P_{hj}(s,t)=\text{Prob}(X(t)=j|X(s)=h),$ (1)

for every possible combination of states $h,j\in S$; $s,t\in[0,\tau),s\leq t$,
where the probability does not depend on the history up to (but not including)
time $s$.

From this we can define transition hazards (or transition intensities) as the
derivative of $P_{hj}(s,t)$ with respect to $t$ (which we assume to exist),
evaluated at $s=t$

$\lambda_{hj}(t)=\lim_{\Delta t\rightarrow 0}\frac{P_{hj}(t,t+\Delta
t)}{\Delta t};$

the cumulative hazard is then defined as

$\Lambda_{hj}(t)=\int\limits_{0}^{t}\\!\lambda_{hj}(u)\,\mathrm{d}u.$

Apart from the definition in formula (1), transition probabilities can also be
defined in a matrix form

$\mathbf{P}(s,t)=\Prodi_{(s,t]}(\mathbf{I}+\boldsymbol{\lambda}(u)\mathrm{d}u).$
(2)

where $\mathbf{P}(s,t)$ represents a square matrix of size $K\times K$. In
such a matrix, every entry represents the transition probability of being in a
given state at time $t$ conditional on being in a given state at time $s$.
Similarly, $\boldsymbol{\lambda}(t)$ is a matrix of size $K\times K$
containing the hazards $\lambda_{hj}(t)$ for all possible transitions
$h\rightarrow j,h\neq j$; whereas diagonal values are defined as
$\lambda_{hh}(t)=-\sum\limits_{j\neq h}\lambda_{hj}(t)$.

Figure 1: Transition diagram for the basic multi-state model.

Both transition hazards and probabilities are of key interest for our work. In
Figure 1 we introduce the multi-state model which will be used as the leading
example in the article. In this model, all patients with a certain disease are
alive at time zero and have no change of their underlying disease status (e.g.
relapse). From the starting state they can then experience two competing
risks: relapse (an intermediate event) or death (i.e. non-relapse mortality,
NRM; an absorbing state). Death can also occur after relapse. We call this
state death after relapse (DaR) which is also an absorbing state. Thus, the
model in Figure 1 is the usual illness-death model where the death state is
split into two different states in order to calculate separate probabilities
of NRM and DaR (which sum up to total mortality). Although the intermediate
event considered is relapse, any other relevant event could be considered,
e.g. recurrence, progression, adverse event, additional treatment, or
recovery. In the following, we assume that the censoring times for relapse and
death are equal, that all relapses before the individual censoring times are
recorded in the dataset and that the event times are observed exactly.

This model contains the key features of a multi-state model: combinations of
competing risks and series of events. All theory and the implementation in
software discussed below are also valid when the model is extended with
additional states. They also hold when individuals start in different states
or are only observed after a certain timepoint (left-truncated observations).

We now turn to estimation. The Nelson-Aalen estimator is the common non-
parametric estimator used for the cumulative transition hazard:

$\widehat{\Lambda}_{hj}(t)=\int\limits_{0}^{t}\\!\frac{\mathrm{d}N_{hj}(u)}{Y_{h}(u)},$
(3)

where $N_{hj}(u)$ represents the counting process that gives the number of
$h\rightarrow j$ transitions in the time interval $[0,u]$, whereas
$Y_{h}(u)=\sum\limits_{i=1}^{n}Y_{h,i}(u)$ gives the size of the at-risk set,
i.e. the number of individuals present in state $h$ just before time $u$
($Y_{h,i}(u)$ indicates whether the ith individual is present in state $h$
just before time $u$).

For estimating the variance of the Nelson-Aalen estimator the Greenwood
estimator is used, as commonly suggested in the literature [16, 17, 18, 19]:

$\widehat{\text{Var}}_{G}(\widehat{\Lambda}_{hj}(t))=\int\limits_{0}^{t}\\!\frac{\mathrm{d}N_{hj}(u)}{Y_{h}(u)(Y_{h}(u)-\mathrm{d}N_{hj}(u))}.$
(4)

For estimating transition probabilities the Aalen-Johansen estimator is
standardly used; it is obtained from formula (2) by replacing
$\boldsymbol{\lambda}(t)$ with the corresponding estimator
$\boldsymbol{\hat{\lambda}}(t)$:

$\mathbf{\hat{P}}(s,t)=\Prodi_{(s,t]}(\mathbf{I}+\boldsymbol{\hat{\lambda}}(u)\mathrm{d}u),$
(5)

or written in discrete form:
$\mathbf{\hat{P}}(s,t)=\displaystyle\prod_{u\in(s,t]}(\mathbf{I}+\Delta\boldsymbol{\widehat{\Lambda}}(u))$,
where $u$ indicates all event times in $(s,t]$ in the dataset.

We point out that both the Nelson-Aalen and Aalen-Johansen estimator are
consistent when the multi-state model is Markov. For a non-Markov process the
Nelson-Aalen estimator is also consistent [20]. We also note that by using the
landmark Aalen-Johansen estimator, each $P_{hj}(s,t)$ can be estimated
consistently [21].

For estimating the variance of transition probabilities we can again use the
Greenwood estimator. We denote it as
$\widehat{\text{Var}}_{G}(\mathbf{\hat{P}}(s,t))$, for which an exact
definition is available in [22, p. 294] and [16].

Here we have provided only the basic tools from multi-state models that are
needed for this work; a more thorough and theoretical introduction can be
found in [23, 24, 25].

###  2.2 Relative survival

We will now introduce the main ideas of relative survival which are relevant
for this work, using the basic survival setting where only one event (death)
is reported in the data. Suppose we are interested in distinguishing between
two causes of death – death due to the disease and its treatment and due to
other (population) causes. To make this distinction we use population
mortality tables as an additional source of information. Such mortality tables
give hazards for every individual based on demographic covariates which have
to be present both in the dataset and in the mortality tables (e.g. age, sex,
year, country) to allow for proper matching. This enables us to evaluate the
amount of deaths due to other (population) causes in the diseased population.

This information has to be appropriately incorporated in a model. We assume
that the (observed) hazard of dying $\lambda(t)$ is the sum of two
(unobserved) hazards, the hazard of dying because of the disease (i.e. the
excess hazard, denoted as $\lambda_{E}(t)$) and the hazard of dying because of
other causes (the population hazard, denoted as $\lambda_{P}(t)$). This
assumption characterizes the additive model:

$\lambda(t)=\lambda_{E}(t)+\lambda_{P}(t).$ (6)

The model is illustrated in Figure 2. In the interest of clarity we give an
exact definition of these hazards. Let $T_{E}$ be the time to excess death and
$T_{P}$ time to population death (i.e. death from other causes). We denote the
minimum of these two times as $T$, the time of death. Time $T$ can be also
subject to a censoring time $T_{C}$ which we assume to be non-informative.
Hence, $\widetilde{T}=\text{min}(T,T_{C})$ is observed, together with the
event indicator at $\widetilde{T}$ (event or censoring). Then the hazard
definition follows

$\lambda(t)=\lim_{\Delta t\to 0}\frac{P(t<T\leq t+\Delta t|T>t)}{\Delta t};$
$\lambda_{E}(t)=\lim_{\Delta t\to 0}\frac{P(t<T_{E}\leq t+\Delta
t|T>t)}{\Delta t};\hskip 15.0pt\lambda_{P}(t)=\lim_{\Delta t\to
0}\frac{P(t<T_{P}\leq t+\Delta t|T>t)}{\Delta t}.$

The two hazards $\lambda_{E}(t)$ and $\lambda_{P}(t)$ are referred to as
cause-specific hazards [26] in analogy to the standard competing risks model
where both causes of death are observed.

It is important to note that in the definition of $\lambda_{E}(t)$ and
$\lambda_{P}(t)$ we condition on $\\{T>t\\}$. This means that we condition on
the observed time and as put in [27] we “stick to the real world”. Another
possibility would be to condition on $\\{T_{E}>t\\}$ and $\\{T_{P}>t\\}$,
respectively, which is commonly done when defining net survival [5]. A more
in-depth study of hazard definitions is provided in [10, Section 7], [5, 9].
We further assume that the population of interest is small compared to the
general population, i.e., that removing the subpopulation of interest from
population life tables would have negligible effect on $\lambda_{P}$.

Figure 2: Illustration of the additive model given in formula (6).

$\Lambda(t)=\int_{0}^{t}\\!\lambda(u)\mathrm{d}u$ can be estimated using the
Nelson-Aalen estimator. $\lambda_{P}(t)$ is estimated using mortality tables:
for the ith individual we can obtain $\lambda_{Pi}(t)$ based on demographic
covariate values. An estimator for the cumulative version of $\lambda_{P}(t)$
is then defined as

$\widehat{\Lambda}_{P}(t)=\int\limits_{0}^{t}\\!\frac{\sum\limits_{i=1}^{n}Y_{i}(u)\mathrm{d}\Lambda_{Pi}(u)}{Y(u)},$
(7)

based on $n$ individuals in the dataset. The definition of
$\widehat{\Lambda}_{P}(t)$ corresponds to the definition of $\Lambda_{P}(t)$
when written with respect to the covariate distribution, see formula (1) in
the Supplementary material where we have also provided a derivation.

The estimator for the cumulative excess hazard $\Lambda_{E}(t)$ is defined as
the difference of the estimators for the observed and population hazard, that
is,

$\widehat{\Lambda}_{E}(t)=\int\limits_{0}^{t}\\!\frac{\mathrm{d}N(u)}{Y(u)}-\int\limits_{0}^{t}\\!\frac{\sum\limits_{i=1}^{n}Y_{i}(u)\mathrm{d}\Lambda_{Pi}(u)}{Y(u)},$
(8)

where $N(u)$ represents the number of events in $[0,u]$.

###  2.3 Extended multi-state model integrating relative survival

####  2.3.1 Estimation of hazards and probabilities

We will now extend this basic multi-state model presented in Section 2.1 with
the relative survival approach introduced in Section 2.2. Thus, we will split
death states in a multi-state model into excess and population counterparts.
To better illustrate the idea, we will extend the example given in Figure 1.
By matching the observed data with population mortality data, we would like to
estimate the amount of population and excess deaths, both before and after the
intermediate event. Such a model is illustrated in Figure 3 – although the
exact cause of death is not present in the dataset (the observed data are the
same as before), we show two additonal arrows and two additional death states
by which we differentiate between excess and population death.

Figure 3: Transition diagram for the extended multi-state model.

We now discuss the quantities of interest and their estimation under this
model. We start with the hazards. For transitions that do not reach death (for
example the transition from Alive relapse-free to Relapse, see Figure 3), the
hazard is the same as in the standard multi-state model (Figure 1). It can
again be estimated by the Nelson-Aalen estimator (3). On the other hand, each
transition reaching death is now split into two (compare transitions 2 and 3
in Figure 1 to the split transitions in Figure 3). The hazards for these
transitions have the same form as in Section 2.2. The estimators follow from
formulas (7), (8) which we now write in general form. For a transition
reaching death $h\rightarrow j$ the estimators of the cumulative hazards for
the population and excess transitions, denoted by $(hj,P)$ and $(hj,E)$, are

$\widehat{\Lambda}_{hj,P}(t)=\int\limits_{0}^{t}\\!\frac{\sum\limits_{i=1}^{n}Y_{h,i}(u)\mathrm{d}\Lambda_{Pi}(u)}{Y_{h}(u)};\widehat{\Lambda}_{hj,E}(t)=\int\limits_{0}^{t}\\!\frac{\mathrm{d}N_{hj}(u)}{Y_{h}(u)}-\int\limits_{0}^{t}\\!\frac{\sum\limits_{i=1}^{n}Y_{h,i}(u)\mathrm{d}\Lambda_{Pi}(u)}{Y_{h}(u)}.$
(9)

The estimators provided in formula (9) have the same form as in formulas (7),
(8), the difference being that now through the at-risk process $Y_{h}(t)$ and
counting process $N_{hj}(t)$ we can take into account additional competing
risks (e.g. relapse) and left-truncation. Transitions from the intermediate
state to death are usually left-truncated since individuals experience the
transition to the intermediate event at different timepoints after time 0.
Partly as a consequence of the later entry in the state and partly because
some demographic features might be associated with an increased risk for the
intermediate event, demographic covariate values of those who experienced the
intermediate event tend to be different from the whole group at baseline.
These aspects are properly taken into account by continuously updating the
population hazard based on the characteristics of the individuals observed to
be at risk for each transition at each timepoint.

Compared to the Nelson-Aalen estimator
$\widehat{\Lambda}_{hj}(t)=\int_{0}^{t}\\!\frac{\mathrm{d}N_{hj}(u)}{Y_{h}(u)}$
which is a step function (jumps occur only at event times),
$\widehat{\Lambda}_{hj,P}(t)$ as defined in (9) is continuous and piecewise
linear, since $\lambda_{Pi}(u)$ is taken from mortality tables. Consequently,
$\widehat{\Lambda}_{hj,E}(t)$ is defined as a difference of a step and a
continuous function (and is càdlàg, like $\widehat{\Lambda}_{hj}(t)$). In
practice, the estimation is done on small intervals (usually on a daily
basis). However, $\widehat{\Lambda}_{hj,P}(t)$ and
$\widehat{\Lambda}_{hj,E}(t)$ are then evaluated at event times only (like the
Nelson-Aalen estimator $\widehat{\Lambda}_{hj}(t)$) and they are taken to be
constant between event times. This procedure causes some bias, which is mostly
negligible. The same approach is taken in [28] where it is also further
discussed (see Section 3.4). Nonetheless, our R implementation allows for
additional evaluation of the estimators at any choice of timepoints (e.g. on a
daily basis). An illustration of the estimation process is given in the
Supplementary material, Section S2.

Having defined the estimates for transition hazards, we now turn to transition
probabilities. Estimator $\mathbf{\hat{P}}(s,t)$ is again defined in matrix
form (5), whereby $\boldsymbol{\hat{\lambda}}(t)$ is extended using excess and
population hazards for all split transitions. The product integral (5) then
provides estimates, where for every split transition $h\rightarrow j$ we
obtain estimates $\widehat{P}_{hj,P}(s,t)$ and $\widehat{P}_{hj,E}(s,t)$,
whereas for all remaining transitions the definition of the estimator of the
transition probabilities stays the same. In practice, transition probabilities
are also estimated on small intervals but evaluated at event times only - a
corresponding example is provided in the Supplementary material.

There are also other measures that are common in relative survival. One such
measure is net survival [5]. However, it is questionable whether this measure
is of interest for a multi-state model as its interpretation becomes less
straightforward.

####  2.3.2 Variance estimation

We consider two approaches for estimating the variances of hazard and
probability estimators - one based on the Greenwood estimator and the second
one based on non-parametric bootstrap (simple resampling with replacement of
individuals).

We start with transition hazards: for non-split transitions we use the usual
Greenwood estimator (4) or bootstrap. To estimate the variances of the split
hazards $\widehat{\Lambda}_{hj,E}(t)$ and $\widehat{\Lambda}_{hj,P}(t)$ we
also consider two options. In the first, we assume that
$\text{Var}(\widehat{\Lambda}_{hj}(t))=\text{Var}(\widehat{\Lambda}_{hj,E}(t))$,
thus we take $\widehat{\text{Var}}_{G}(\widehat{\Lambda}_{hj}(t))$ as the
estimator of $\text{Var}(\widehat{\Lambda}_{hj,E}(t))$. In this case, one
assumes that the population hazard $\widehat{\Lambda}_{hj,P}(t)$ is fixed and
has variance zero. This assumption is often made because of the use of
mortality tables which are deterministic. However, as is clear from formula
(9), the $\widehat{\Lambda}_{hj,P}(t)\text{'s}$ are not deterministic since
$Y_{h}(t)$ brings some variability through the composition of the study
population in terms of the demographic covariates (age, sex, year etc.). As a
second option for evaluating variability we use the bootstrap, which will help
investigate if this additional assumption is reasonable.

Analogously, the variance of estimates for the transition probabilities is
estimated using the Greenwood and bootstrap options.

####  2.3.3 Confidence interval methods

We consider four confidence interval (CI) methods for transition hazards and
probabilities. Although the model of interest is still the model introduced in
the previous subsection, the methods are general and can be applied in any
non-parametric multi-state Markov model. We denote by $\theta$ the transition
hazard or probability at a given time for a given transition and
$\hat{\theta}$ its corresponding estimate. By denoting
$\widehat{\text{Var}}_{G}(\hat{\theta})$ and
$\widehat{\text{Var}}_{boot}(\hat{\theta})$ to be the Greenwood and bootstrap
estimators of the variance, respectively, we define the four CI methods.

(CI1)

Plain scale using Greenwood (plain.G): a symmetrical $(1-\alpha)\%$ CI on the
plain scale using the Greenwood estimator for the variance:

$\hat{\theta}\pm\Phi(1-\frac{\alpha}{2})\cdot\sqrt{\widehat{\text{Var}}_{G}(\hat{\theta})},$
(10)

where $\Phi$ denotes the cdf of the standard normal distribution.

(CI2)

Plain scale using bootstrap (plain.boot): a symmetrical $(1-\alpha)\%$ CI on
the plain scale using the bootstrap estimator for the variance:

$\hat{\theta}\pm\Phi(1-\frac{\alpha}{2})\cdot\sqrt{\widehat{\text{Var}}_{boot}(\hat{\theta})}.$
(11)

(CI3)

Log scale using bootstrap (log.boot): a $(1-\alpha)\%$ CI on the log scale
(based on the delta method) using the bootstrap estimator for the variance:

$\hat{\theta}\cdot\text{exp}\Big{\\{}\pm\Phi(1-\frac{\alpha}{2})\cdot\frac{\sqrt{\widehat{\text{Var}}_{boot}(\hat{\theta})}}{\hat{\theta}}\Big{\\}}.$
(12)

(CI4)

Quantiles using bootstrap (q.boot): a $(1-\alpha)\%$ CI where the limits of
the CI are taken to be the $2.5\%$ and $97.5\%$ quantiles of $\hat{\theta}$
over the bootstrap samples.

##  3 Simulations

We analyse the characteristics of the developed methodology presented in
Section 2.3 through simulations. First, in Section 3.1 we introduce all
details regarding how the simulations were performed and in Section 3.2 we
present the results. All simulations were done in R, version 3.5.2 [29] using
the packages mstate [30] and relsurv [31]. The main simulation functions have
been added in a separate .R-file (simulation_functions.R) in the Supplementary
material.

###  3.1 Defining the simulation study

In this section we present the aims, data-generating mechanisms, estimands,
methods and performance measures (i.e. ADEMP) used for simulation. These are
reported as proposed in the tutorial by Morris et al. [32].

####  3.1.1 Aims

Our main goal is to investigate how the newly proposed estimators for the
extension of multi-state models with population mortality perform in a wide
range of simulation scenarios. We will evaluate their bias, standard errors
(SE) and coverage probabilities.

####  3.1.2 Data-generating mechanisms

Transition diagram: we use the same transition diagram as in Figure 3. At time
$0$, all individuals are alive event-free. All individuals are followed up to
a maximum of $10$ years.

Generating event times: for every individual we generate a trajectory, where
the variables and event times are obtained in the following way:

  1. 1.

Demographic variables $D$ are simulated independently for every individual:
age (uniform), sex (binary), year (uniform). In all simulations, sex is
distributed Bernoulli($0.5$) and year of diagnosis is taken
Uniform[1990-01-01, 2000-01-01]. Age is also simulated with the Uniform
distribution on a prespecified interval that will be adjusted for every
simulating scenario so that the amount of population deaths is controlled.

  2. 2.

We first have three event times for the following competing risks: relapse,
population NRM and excess NRM. In practice, we only observe NRM and relapse,
but in the simulations we also generate the underlying population and excess
times to death. If an individual reaches state Relapse, times to population
DaR and excess DaR are then simulated conditional on the time of relapse. All
events times are simulated independently [33].

  3. 3.

Times to death from population causes are simulated using population mortality
tables for Slovenia (which are available in R \- package relsurv, object
slopop). For simulating the transition from Relapse to population DaR the
demographic covariate values at time of relapse are considered.

  4. 4.

Remaining event times (time to relapse, excess NRM and excess DaR) are
simulated using
$\lambda_{i}(t)=\lambda_{0}(t)\cdot\exp(\boldsymbol{\beta}^{\top}D_{i})$.

     * •

In the simplest case, we take the hazard to be constant (i.e. we use the
exponential distribution). In this case, the baseline hazard is constant and
there are no covariate effects ($\lambda_{0}(t)=\lambda_{0}$).

     * •

In some simulating scenarios, we take Weibull distributed event times for
these transitions. In this case, $\lambda_{0}(t)$ is the hazard of a Weibull
distribution with parameters $a$ (rate) and $b$ (shape);
$\lambda_{0}(t)=abt^{b-1}$. In this case, no covariate effects are considered.
We further note that event times for the transition from relapse to excess DaR
are simulated using a left-truncated Weibull distribution [34] where one
conditions on the time of relapse; by this the Markov assumption is fulfilled.

     * •

Furthermore, we take simulation scenarios with a covariate effect (where
$\boldsymbol{\beta}$ is not equal to zero). The covariate effect is included
with respect to age, where age is first centered around 0.

  5. 5.

Although we simulate event times for both population and excess transitions,
only the minimum of these is included in the dataset as observed value.

  6. 6.

The generated times are then subjected to censoring. Censoring times are
simulated from the exponential distribution. The rate parameter of the
distribution is modified for every scenario such that approximately 20% of all
individuals are censored at the end of follow-up time ($10$ years).

Simulation scenarios: as mentioned, all population-related event times are
generated using the same mortality tables, but population hazards vary between
scenarios by adjusting the age distribution. Based on how the remaining event
times are distributed and the magnitude of their hazards we define the five
following simulation scenarios:

  1. 1.

$exp.small$: for all transitions that are not population-related, the
exponential distribution is used and the amount of population deaths is small.

  2. 2.

$exp.large$: for all transitions that are not population-related, the
exponential distribution is used and the amount of population deaths is close
to the amount of excess deaths.

  3. 3.

$weibull$: for all transitions that are not population-related, the Weibull
distribution is used and the amount of population deaths is small.

  4. 4.

$cov.eff.pos$: we allow for a covariate effect for certain transitions, the
baseline hazard is again constant. The covariate effects of age on the hazard
of having excess NRM and relapse are positive. No covariate effect is taken
for the hazard of having excess DaR.

  5. 5.

$cov.eff.neg$: we allow for a covariate effect for certain transitions, the
baseline hazard is again constant. The covariate effect of age on the hazard
of having excess NRM is positive, whereas on relapse it is negative. No
covariate effect is taken for the hazard of having excess DaR.

Details on exact parameters can be found in the Supplementary material. We
additionally assess the impact of sample size on the estimators by varying
sample size: $n=500,1000$ and $2000$ subjects.

####  3.1.3 Estimands

Our main measures of interest are the transition hazards and probabilities of
the multi-state model. All values are assessed at $1,2,5$ and $10$ years,
together with the corresponding standard errors and confidence intervals. All
these measures have been defined in Section 2.3.

Exact transition hazards and probabilities (the true values) were calculated
by integrating with respect to the covariate distribution $D$, which was not
straightforward since times to population death were generated using mortality
tables (see Section S1 in the Supplementary material). The numerical
integration was done in R.

####  3.1.4 Methods

We use the non-parametric approach outlined in Section 2.3 to estimate all
transition hazards and probabilities, and we compare the performance of the
variance estimators and the four different confidence interval methods.

####  3.1.5 Performance measures

Bias: let $\theta$ be the transition hazard or probability estimated at a
given time for a given transition, whereas $\hat{\theta}_{i}$ is the estimate
of $\theta$ for the ith replication of a simulation scenario. We define

$\overline{\hat{\theta}}=\frac{1}{n_{sim}}\sum_{i=1}^{n_{sim}}\hat{\theta}_{i}.$

Then the definition of absolute and relative bias follows:

$\text{absolute bias}=\theta-\overline{\hat{\theta}},\hskip
60.0pt\text{relative bias}=\frac{\theta-\overline{\hat{\theta}}}{\theta}.$

Relative bias is reported in the simulation results. This allows for an easy
comparison of results across different timepoints.

Standard errors: we compare the two approaches for estimating variances
(Greenwood and bootstrap) as defined in Section 2.3.2. For the bootstrap $100$
replications are taken. We also calculate empirical standard errors, i.e. the
standard deviation across all simulation replications:

$\widehat{\textrm{SE}}=\sqrt{\frac{1}{n_{\text{sim}}-1}\sum_{i=1}^{n_{\text{sim}}}(\hat{\theta}_{i}-\overline{\hat{\theta}})^{2}}.$

Coverage probabilities: $95\%$ confidence intervals (CI) will be calculated
for transition hazards and probabilities using the four methods introduced in
Section 2.3.3. The estimated coverage probability (CP) is then the proportion
of times the resulting CI contains the true estimand value.

We also use the coverage probability for calculating the number of
replications $n_{\text{sim}}$ needed in the simulations. As in Section 5.3 in
Morris et al. [32], we require that the SE for the coverage probability equals
$0.5\%$. Knowing that the expected coverage probability equals $95\%$, it
follows that $n_{\text{sim}}=1900$. We note that if we would require the SE
for the coverage probability to be $1\%$ we would obtain $n_{sim}=475$. In our
simulations we round this to $n_{\text{sim}}=2000$.

###  3.2 Results

####  3.2.1 Bias

Bias was negligible in all scenarios. Figure 4 shows the relative bias for
cumulative transition hazards for all scenarios evaluated at the end of
follow-up ($10$ years). An equivalent figure for transition probabilities is
given in Figure S4 in the Supplementary material. On average, relative bias
was smaller when the sample size increased, which is expected. It seems that
the relative bias was distributed around 0 and that behaviour was appropriate
for all transitions and scenarios. For simulation scenarios with a covariate
effect a slightly larger deviation from $0$ was seen, but when the sample size
was large enough (e.g. $2000$), the relative bias was close to $0$.

Figure 4: Relative bias for cumulative hazards evaluated at end of follow-up
($10$ years). States: Alive relapse-free (ARF), Relapse, population NRM
(NRM.p), excess NRM (NRM.e), population DaR (DaR.p), excess DaR (DaR.e).

####  3.2.2 Standard errors

In Figure 5 we show standard errors for scenario $exp.large$ with sample size
$2000$. Through this example we will explain the crucial characteristics that
were present in all of the simulation scenarios. The Greenwood option gave a
zero estimate of the variance for population-related hazards and consequently,
the estimated variance was also smaller for the corresponding transition
probabilities. The bootstrap option better reflected the uncertainty in the
estimates, as is confirmed by its values being close to the empirical SE.
Furthermore, standard errors were smaller for population-related transitions
than for their excess counterparts. When the amount of population deaths was
increased, the SE of these population-related transitions increased as well.

Figure 5: Standard erorrs for cumulative hazards and transition probabilities
for scenario $exp.large$ with sample size $2000$ (average over 2000
replications). States: Alive relapse-free (ARF), Relapse, population NRM
(NRM.p), excess NRM (NRM.e), population DaR (DaR.p), excess DaR (DaR.e).

When looking at the left-truncated transition going from relapse to excess
DaR, we can see that both variance estimators for the cumulative hazards,
Greenwood and bootstrap-based, gave evidently smaller SEs than the empirical
SE. This occurred for all simulation scenarios. As further explained in
Section S3 of the Supplementary material, this phenomenon does not only appear
for the extended multi-state model but is already present in a basic multi-
state model for left-truncated transitions. Solving this problem is out of the
scope of this paper.

For transition probabilities (Figure 5) both estimator options were close to
the empirical SE, the population-related transitions being an exception. In
this case, the bootstrap was again a reliable option. The standard errors for
population-related transition probabilities increased when the amount of
population deaths increased.

We find the bootstrap-based method a reliable option for estimating variances
for all transitions in the wide range of scenarios we have explored. For some
simulation scenarios (those with a covariate effect) there was a small
difference between the bootstrap estimate and the empirical SE for samples of
size $500$ and $1000$. However, these differences diminished for samples of
size $2000$.

We additionally compared the standard errors for the split transitions to the
standard errors for the observed transitions. In most cases, the empirical
standard errors for the observed transitions were very similar to the ones for
the excess-related transitions. This has been proven to hold asymptotically
for transition hazards, i.e.
$\text{Var}(\widehat{\Lambda}_{hj}(t))\approx\text{Var}(\widehat{\Lambda}_{hj,E}(t))$
[22, p. 202]. On the other hand, standard errors for population-related
transitions were smaller than for their excess-counterparts and they were
mostly dependent upon the amount of population deaths and variability in the
distribution of the demographic variables.

####  3.2.3 Confidence intervals

Figure 6 shows coverage probabilities for transition probabilities at the end
of follow-up. An equivalent graph for transition hazards is shown in the
Supplementary material (Figure S5). On average, coverage probabilities were
more stable and closer to the nominal value for transition probabilities than
for hazards. By increasing the sample size the coverage probabilities also
tended to increase and stay close to the nominal value (results not shown).

Figure 6: Coverage probabilities for transition probabilities evaluated at the
end of follow-up ($10$ years) (averaged across sample sizes). Method $plain.G$
gives inadequate coverage probabilities for population NRM and population DaR,
many of which are smaller than $85\%$ and thus are not shown in the figure.
States: Alive relapse-free (ARF), Relapse, population NRM (NRM.p), excess NRM
(NRM.e), population DaR (DaR.p), excess DaR (DaR.e).

In almost all cases, $q.boot$ gave smaller CPs than the other methods. For
population-related transitions, the $plain.G$ method gave inadequate
confidence intervals - on the hazard level we got zero-width CIs, on the
probability level they were somewhat larger but still far from the nominal
value $95\%$ (many of them were smaller than $85\%$ and thus not shown in
Figure 6). This is a consequence of the variance estimator that was taken to
equal $0$ for these transitions.

The remaining two methods $plain.boot$ and $log.boot$ which were based on
bootstrap SE estimates behaved similarily in the simuations, but the latter
method gave much better results for some cases. For left-truncated and
population-related transitions the $log.boot$ method gave larger CPs which
were closer to the nominal coverage - this was especially evident at earlier
timepoints. An example of this is shown in Figure S6 in the Supplementary
material. Thus, $log.boot$ is the preferred choice between the four introduced
CI methods.

However, there were cases where even $log.boot$ gave lower coverage
probabilities than expected; this was most evident for some left-truncated and
population-related transitions for scenarios $cov.eff.pos$ and $cov.eff.neg$.
This phenomenon occurred whenever the bootstrap SE estimate was lower than the
empirical SE.

##  4 Illustration

The new methodology is briefly illustrated by a reanalysis of the data
presented in [4], which is taken as a motivating example in this article. The
study population consisted of patients who had received a first allogeneic
hematopoietic stem cell transplantation (alloHCT) for MDS or sAML between
January 2000 and December 2012 and were recorded in the registry of the EBMT
[35]. Further details on the selection criteria, patient characteristics and
outcomes are available in [4].

Patients were followed since alloHCT and possible outcomes were
relapse/progression (for the remainder of this section we call this relapse)
or death where we again distinguish between non-relapse mortality (NRM) and
death after relapse (DaR). Thus, a multi-state model with a transition diagram
as in Figure 1 is applicable, which can be further extended as in Figure 3. To
better illustrate the usefulness of the methodology proposed, a subsample from
the original dataset was taken, only including patients aged $\geq 60$ years
at alloHCT and still alive relapse-free at 2 years after alloHCT (the landmark
time). Population mortality is especially relevant for this group of older
patients who have survived the first hazardous period after alloHCT. Data from
a total of $753$ patients from $19$ European countries were included, with a
median age of $63.9$ years, where $62\%$ of them were male and $38\%$ female.
For this example, we restrict follow-up time to $8$ years.

Figure 7 shows the estimated cumulative hazards (left graph) and transition
probabilities (right graph) for the extended multi-state model, whereas in the
Supplementary material, these estimated values are shown with the
corresponding confidence intervals, see Figure S7 and S8. The huge cumulative
hazard of excess DaR is a sign of the very poor prognosis after relapse. The
cumulative hazard of experiencing relapse is for most of the time slightly
larger than the cumulative hazard of excess NRM. The hazards for the
population-related transitions are very similar, their difference only caused
by the different composition of the risk sets in terms of demographic
variables in the states Alive relapse-free and Relapse.

Figure 7: Estimated cumulative hazards (left) and transition probabilities
(right) for patients 65 years or older at alloHCT and alive relapse-free at
the 2-year landmark time estimated in the EBMT dataset (as introduced in
Section 4). Estimates provided up to $8$ years since the landmark time.
States: Alive relapse-free (ARF), Relapse, population NRM (NRM.p), excess NRM
(NRM.e), population DaR (DaR.p), excess DaR (DaR.e).

The transition probabilities (right graph in Figure 7) provide additional
insight into the outcomes of the transplanted patients. The probability of
being alive and relapse-free steadily decreases through time and reaches a
probability of $56.9\%$ ($95\%$ $log.boot$ CI: $[52.1,62.1]$) at $8$ years
(all probabilities in this paragraph are reported at this time). The pattern
shown by the hazards is propagated by the transition probabilities. The
probability of dying due to population causes (population NRM or population
DaR) equals $10.9\%$. The probability of population NRM equals $10.3\%$
($95\%$ CI: $[9.7,10.9]$), whereas its excess counterpart (the probability of
excess NRM) is $15.1\%$ ($95\%$ CI: $[11.3,20.2]$). On the other hand, the
probability of population DaR is much smaller than of excess DaR: $0.6$
($95\%$ CI: $[0.4,1.0]$) and $13.3$ ($95\%$ CI: $[10.7,16.5]$), respectively,
which is due to the strong competing risk of excess DaR. The probability of
being in the relapse state (and having no further event) is steady after two
years reaching $3.9\%$ ($95\%$ CI: $[1.7,8.8]$) at end of follow-up.

The clinical interpretation of these results is that death after relapse is
almost identical to death due to relapse. In the past, only young patients
were referred for alloHCT. For these younger patients (below the age of 60
years at the time of alloHCT) non-relapse mortality is almost exclusively
treatment-related mortality (where treatment refers both to the alloHCT and to
other treatments before and after). Yet age at HCT has continuously increased
in the past decades, and the presumption that non-relapse mortality equals
treatment-related mortality is not valid for the older transplanted
population. In our example, this is demonstrated by the fact that $40.6\%$ of
NRM at 10 years after alloHCT for this landmark cohort can be attributed to
population mortality. Thus, their treatment-related mortality can better be
estimated by excess NRM.

The R package mstate [30] allows for the introduced extension. After the
transition hazards are estimated using the function msfit, they are supplied
to the newly developed msfit.relsurv function. This produces a new and
extended msfit object containing the hazards for the excess and population
transitions and their variances, calculated by the methods described in
Section 2.3.2. The extended msfit object can then be supplied to function
probtrans, through which transition probability estimates and their variance
estimates are obtained. Furthermore, plain and log confidence intervals (as
defined in Section 2.3.3) can be obtained using the summary.msfit and
summary.probtrans functions. Thus the original package has been upgraded in
such a way that only a small additional step is needed when coding (function
msfit.relsurv). A full code example using mstate is provided in Section S4 of
the Supplementary material.

##  5 Discussion

In this paper, we have introduced an extension of a multi-state model
containing death states where the exact cause of death is not known. We were
able to account for death by other (population) causes in such a model using
external population mortality tables. The relative survival theory has allowed
for such an extension where we have solely focused on non-parametric
estimation in a Markov setting. Thus, we have successfully integrated the
basic relative survival idea (of splitting a death state in excess and
population counterparts) in a general multi-state framework for any transition
reaching death.

The present work focuses on transition hazards (through which the extension
has been done) and transition probabilities, the latter being more easily
interpretable for a non-statistical audience. The proposed estimators are
defined in such a way that the information from mortality tables is properly
taken into account. Based on our simulation study we can conclude that as
expected, estimation bias is very minor.

By using the extended multi-state model we also have to assume the additive
model (formula (6)) for all split transitions. Thus, in order for the excess
hazard to be a proper hazard, it needs to be positive across the whole time
interval. There are many practical settings where this does not hold (e.g.
[36]) because the population under study has a better life expectancy than the
general population, e.g., because primarily patients with a more robust
general health might survive certain diagnoses. If the model is still applied
even though it does not really hold, negative cumulative hazards and
transition probabilities are obtained, which could still give a useful
message, if interpreted with caution. In such a case, negative cumulative
hazards can only be obtained for excess-related transitions which then affect
the transition probabilities of the excess and population death states.

Throughout the paper, we only split transitions to death states. However, the
methods (together with the software implementation) allow transitions reaching
intermediate events to be split as well. In this case, we would need adequate
population tables for those specific transitions. For example, if we want to
examine if patients with an autoimmune disease run more risk of a Covid-19
infection and death thereafter than the general population, a model could be
developed with starting state Diagnosis, intermediate event Covid-19 infection
and absorbing states Death before and after Covid-19 infection. Based on
observed data for a cohort of these patients and population tables for the
hazard of infection by Covid-19 and death after infection, the transition to
the intermediate event could be split into a population and an excess part and
similarly the transition from Infection to Death after infection could be
split. A different example is available in [15] where a model for Hodgkin
lymphoma patients with a disease of the circulatory system (DCS) as an
intermediate event is considered.

The suggested theory only uses a small part of the field of relative survival;
it does not consider any other relative survival measures, e.g. net survival
or the relative survival ratio, and their extension to a general multi-state
model. We leave this point for further research. We note that the crude
probability of death (a common measure in relative survival) has its analogon
in the split transition probabilities $\widehat{P}_{hj,P}(s,t)$ and
$\widehat{P}_{hj,E}(s,t)$.

For estimating the variance of the proposed measures we propose the bootstrap
approach as best which was a suitable option in most of the simulation
scenarios. The Greenwood estimator was also considered as a closed-form
estimator. This required the additional assumption for the variance estimates
of population-related transitions to be $0$. This assumption proved to be
inadequate since the population-related transitions have substantial
variability, whereas for the remaining transitions the Greenwood option was
suitable in most cases. Instead of making such an assumption, a non-parametric
variance estimator for population-related transitions would be of great use.
However, such an estimator has not been developed (as far as we know) and thus
this additional assumption was required. Developing such a variance estimator
for the population-related transitions could find its use for many measures
related to relative survival.

An additional problem with variance estimation of left-truncated transitions
has been also pointed out (and further illustrated in the Supplementary
material). This problem is inherent to a basic multi-state model already and
is out of the scope of this paper but can be tackled in future work.

In the simulation study, we have also considered four confidence interval
methods. The $log.boot$ performed best, and is therefore our recommended
option. Even though there were cases when the coverage probability of
$log.boot$ was evidently smaller than the nominal value, this method proved to
be a reliable choice across simulations. Nevertheless, there is still room for
improvement, either by defining a new confidence interval method or taking an
even better variance estimator.

All work has been implemented in R and is free for use through the mstate
package [30]. The implementation is based on functions from package relsurv
[31] which has allowed for an easy and optimal usage of mortality tables which
is otherwise not trivial. The current R implementation also allows for non-
parametric estimation based on left-truncated observations in a standard
survival setting (which has not been always considered in previous software
implementations) and for future extensions to Cox-model based semi-parametric
models. From what we know this is the first such implementation in R. Its use
has been illustrated through an application. The Supplementary material shows
self-contained code that can help the reader to apply this methodology.

To conclude, we believe that the proposed extension of a non-parametric multi-
state model could find its use in a wide range of practical applications.
Especially for interventions encompassing potential long-term beneficial and
adverse effects offered to older patients with a remaining life-expectancy of
several years or even a decade, taking into account population mortality as
part of the different components of mortality is relevant. In this paper we
have mostly focused on the theoretical aspects of this work. By conducting a
thorough simulation study, we have further understood the properties of the
estimators and come across additional challenges, some of which have to be
dealt with in the future (e.g. variance estimation for left-truncated
transitions, a theoretical variance estimator for population-related
transitions, an improved confidence interval method). We have devised an R
application that relies on previous multi-state implementations in R. Thus by
providing such a tool we believe that this work can be beneficial in many
practical settings.

### Acknowledgements

Damjan Manevski is a young researcher funded by the Slovene Research Agency
(ARRS). The authors acknowledge the financial support by ARRS (Methodology for
data analysis in medical sciences, P3-0154 and project J3-1761). We thank the
EBMT Data Office in Leiden for providing the MDS dataset.

### Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * Mariotto et al. [2014] Angela B. Mariotto, Anne-Michelle Noone, Nadia Howlader, Hyunsoon Cho, Gretchen E. Keel, Jessica Garshell, Steven Woloshin, and Lisa M. Schwartz.  Cancer Survival: An Overview of Measures, Uses, and Interpretation.  _Journal of the National Cancer Institute. Monographs_ , 2014(49):145–186, November 2014.  ISSN 1052-6773.  doi: 10.1093/jncimonographs/lgu024.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4829054/.
  * Begg and Schrag [2002] Colin B. Begg and Deborah Schrag.  Attribution of deaths following cancer treatment.  _Journal of the National Cancer Institute_ , 94(14):1044–1045, July 2002.  ISSN 0027-8874.  doi: 10.1093/jnci/94.14.1044.
  * Percy et al. [1981] C. Percy, E. Stanek, and L. Gloeckler.  Accuracy of cancer death certificates and its effect on cancer mortality statistics.  _American Journal of Public Health_ , 71(3):242–250, March 1981.  ISSN 0090-0036.  doi: 10.2105/ajph.71.3.242.
  * Schetelig et al. [2019] Johannes Schetelig, Liesbeth C. de Wreede, Michel van Gelder, Linda Koster, Jürgen Finke, Dietger Niederwieser, Dietrich Beelen, G. J. Mufti, Uwe Platzbecker, Arnold Ganser, Silke Heidenreich, Johan Maertens, Gerard Socié, Arne Brecht, Matthias Stelljes, Guido Kobbe, Liisa Volin, Arnon Nagler, Antonin Vitek, Thomas Luft, Per Ljungman, Ibrahim Yakoub-Agha, Marie Robin, and Nicolaus Kröger.  Late treatment-related mortality versus competing causes of death after allogeneic transplantation for myelodysplastic syndromes and secondary acute myeloid leukemia.  _Leukemia_ , 33(3):686–695, March 2019.  ISSN 1476-5551.  doi: 10.1038/s41375-018-0302-y.  URL https://www.nature.com/articles/s41375-018-0302-y.
  * Pohar Perme et al. [2012] M. Pohar Perme, J. Stare, and J. Estève.  On estimation in relative survival.  _Biometrics_ , 68(1):113–120, 2012.
  * HMD [2016] Human mortality database, 2016.  Available at www.mortality.org or www.humanmortality.de.
  * Ederer et al. [1961] F. Ederer, L. M. Axtell, and S. J. Cutler.  _The relative survival rate: a statistical methodology_ , volume 6, pages 101–121.  National Cancer Institute Monograph, 1961.
  * Hakulinen [1982] T. Hakulinen.  Cancer survival corrected for heterogeneity in patient withdrawal.  _Biometrics_ , 38:933–942, 1982.
  * Sasieni and Brentnall [2016] Peter Sasieni and Adam R. Brentnall.  On standardized relative survival.  _Biometrics_ , 73(2):473–482, 2016.  ISSN 1541-0420.  doi: 10.1111/biom.12578.  URL http://dx.doi.org/10.1111/biom.12578.
  * Pavlič and Pohar Perme [2018] Klemen Pavlič and Maja Pohar Perme.  Using pseudo-observations for estimation in relative survival.  _Biostatistics_ , 2018.  doi: 10.1093/biostatistics/kxy008.  URL +http://dx.doi.org/10.1093/biostatistics/kxy008.
  * Huszti et al. [2012] Ella Huszti, Michal Abrahamowicz, Ahmadou Alioum, Christine Binquet, and Catherine Quantin.  Relative survival multistate Markov model.  _Statistics in Medicine_ , 31(3):269–286, February 2012.  ISSN 1097-0258.  doi: 10.1002/sim.4392.
  * Belot et al. [2011] Aurélien Belot, Laurent Remontet, Guy Launoy, Valérie Jooste, and Roch Giorgi.  Competing risk models to estimate the excess mortality and the first recurrent-event hazards.  _BMC Medical Research Methodology_ , 11(1):78, May 2011.  ISSN 1471-2288.  doi: 10.1186/1471-2288-11-78.  URL https://doi.org/10.1186/1471-2288-11-78.
  * Gillaizeau et al. [2017] Florence Gillaizeau, Etienne Dantan, Magali Giral, and Yohann Foucher.  A multistate additive relative survival semi-Markov model.  _Statistical Methods in Medical Research_ , 26(4):1700–1711, August 2017.  ISSN 1477-0334.  doi: 10.1177/0962280215586456.
  * Crowther and Lambert [2017] Michael J. Crowther and Paul C. Lambert.  Parametric multistate survival models: Flexible modelling allowing transition-specific distributions with application to estimating clinically useful measures of effect differences.  _Statistics in Medicine_ , 36(29):4719–4742, 2017\.  ISSN 1097-0258.  doi: https://doi.org/10.1002/sim.7448.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.7448.
  * [15] Caroline E. Weibull, Paul C. Lambert, Sandra Eloranta, Therese M. L. Andersson, Paul W. Dickman, and Michael J. Crowther.  A multistate model incorporating estimation of excess hazards and multiple time scales.  _Statistics in Medicine_ , n/a(n/a).  ISSN 1097-0258.  doi: https://doi.org/10.1002/sim.8894.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.8894.
  * de Wreede et al. [2010] Liesbeth C. de Wreede, Marta Fiocco, and Hein Putter.  The mstate package for estimation and prediction in non- and semi-parametric multi-state and competing risks models.  _Computer Methods and Programs in Biomedicine_ , 99(3):261–274, 2010.  ISSN 0169-2607.  doi: https://doi.org/10.1016/j.cmpb.2010.01.001.  URL https://www.sciencedirect.com/science/article/pii/S0169260710000027.
  * Allignol et al. [2010] Arthur Allignol, Martin Schumacher, and Jan Beyersmann.  A Note on Variance Estimation of the Aalen–Johansen Estimator of the Cumulative Incidence Function in Competing Risks, with a View towards Left-Truncated Data.  _Biometrical Journal_ , 52(1):126–137, 2010.  doi: 10.1002/bimj.200900039.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.200900039.
  * Friedrich et al. [2017] Sarah Friedrich, Jan Beyersmann, Ursula Winterfeld, Martin Schumacher, and Arthur Allignol.  Nonparametric estimation of pregnancy outcome probabilities.  _The Annals of Applied Statistics_ , 11:840–867, 06 2017\.  doi: 10.1214/17-AOAS1020.
  * Rousson et al. [2019] Valentin Rousson, Arthur Allignol, Alexandre Aurousseau, Ursula Winterfeld, and Jan Beyersmann.  Stabilizing cumulative incidence estimation of pregnancy outcome with delayed entries.  _Biometrical Journal_ , 61(5):1290–1302, 2019\.  doi: 10.1002/bimj.201700237.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.201700237.
  * Datta and Satten [2001] S. Datta and G. A. Satten.  Validity of the Aalen-Johansen estimators of stage occupation probabilities and Nelson-Aalen estimators of integrated transition hazards for non–Markov models.  _Statistics and Probability Letters_ , 55:403–411, 2001\.
  * Putter and Spitoni [2018] Hein Putter and Cristian Spitoni.  Non-parametric estimation of transition probabilities in non-Markov multi-state models: The landmark Aalen–Johansen estimator.  _Statistical Methods in Medical Research_ , 27(7):2081–2092, July 2018.  ISSN 0962-2802.  doi: 10.1177/0962280216674497.  URL https://doi.org/10.1177/0962280216674497.
  * Andersen et al. [1993] P. K. Andersen, O. Borgan, R. D. Gill, and N. Keiding.  _Statistical models based on counting processes_.  Springer-Verlag, New York, 1993.
  * Andersen and Keiding [2002] P. K. Andersen and N. Keiding.  Multi-state models for event history analysis.  _Statistical Methods in Medical Research_ , 11:91–115, 2002\.
  * Putter et al. [2007] H. Putter, M. Fiocco, and R. B. Geskus.  Tutorial in biostatistics: competing risks and multi-state models.  _Statistics in Medicine_ , 26:2389–2430, 2007.
  * Klein et al. [2014] John P. Klein, Hans C. van Houwelingen, Joseph G. Ibrahim, and Thomas H. Scheike, editors.  _Handbook of survival analysis_.  Chapman & Hall/CRC handbooks of modern statistical methods. CRC Press, Boca Raton, FL, 2014.  ISBN 9781466555662.
  * Tsiatis [2005] A. A. Tsiatis.  _Competing Risks_.  wiley, 2005.  ISBN 9780470011812.  doi: 10.1002/0470011815.b2a03035.  URL http://dx.doi.org/10.1002/0470011815.b2a03035.
  * Andersen and Keiding [2012] Per Kragh Andersen and Niels Keiding.  Interpretability and importance of functionals in competing risks and multistate models.  _Statistics in Medicine_ , 31(11-12):1074–1088, 2012.  ISSN 1097-0258.  doi: 10.1002/sim.4385.  URL http://dx.doi.org/10.1002/sim.4385.
  * Pohar Perme [2017] M. Pohar Perme.  _relsurv: Relative Survival_ , 2017.  URL http://CRAN.R-project.org/package=relsurv.  R package version 2.1-1.
  * R Core Team [2018] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2018.  URL https://www.R-project.org/.
  * de Wreede et al. [2011] Liesbeth de Wreede, Marta Fiocco, and Hein Putter.  mstate: An R Package for the Analysis of Competing Risks and Multi-State Models.  _Journal of Statistical Software, Articles_ , 38(7):1–30, 2011.  ISSN 1548-7660.  doi: 10.18637/jss.v038.i07.  URL https://www.jstatsoft.org/v038/i07.
  * Perme and Pavlic [2018] Maja Pohar Perme and Klemen Pavlic.  Nonparametric Relative Survival Analysis with the R Package relsurv.  _Journal of Statistical Software, Articles_ , 87(8):1–27, 2018.  ISSN 1548-7660.  doi: 10.18637/jss.v087.i08.  URL https://www.jstatsoft.org/v087/i08.
  * Morris et al. [2019] Tim P. Morris, Ian R. White, and Michael J. Crowther.  Using simulation studies to evaluate statistical methods.  _Statistics in Medicine_ , 38(11):2074–2102, 2019\.  doi: 10.1002/sim.8086.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.8086.
  * Meira-Machado and Faria [2014] L. Meira-Machado and S. Faria.  A Simulation Study Comparing Modeling Approaches in an Illness-Death Multi-State Model.  _Communications in Statistics - Simulation and Computation_ , 43(5):929–946, 2014.  doi: 10.1080/03610918.2012.718841.  URL https://doi.org/10.1080/03610918.2012.718841.
  * Wingo [1989] Dallas R. Wingo.  The left-truncated Weibull distribution: theory and computation.  _Statistical Papers_ , 30(1):39–48, 1989.  doi: 10.1007/bf02924307.
  * [35] The EBMT Patient Registry.  URL https://www.ebmt.org/ebmt-patient-registry.
  * Antero-Jacquemin et al. [2018] Juliana Antero-Jacquemin, Maja Pohar-Perme, Grégoire Rey, Jean-François Toussaint, and Aurélien Latouche.  The heart of the matter: years-saved from cardiovascular and cancer deaths in an elite athlete cohort with over a century of follow-up.  _European Journal of Epidemiology_ , 33(6):531–543, June 2018.  ISSN 1573-7284.  doi: 10.1007/s10654-018-0401-0.