Time-dynamic inference for non-Markov transition probabilities under independent right-censoring

TL;DR

This paper establishes weak convergence results for non-Markov transition probability estimators under right-censoring, enabling more accurate inference and confidence band construction in complex survival analysis scenarios.

Contribution

It proves weak convergence for a general non-Markov estimator and its bootstrap, providing new tools for inference in non-Markov processes with right-censoring.

Findings

Bootstrap methods enable resampling-based inference.

Confidence bands reveal biases in Markov-based estimators.

Application to liver cirrhosis data shows different conclusions from Markov assumptions.

Abstract

In this article, weak convergence of the general non-Markov state transition probability estimator by Titman (2015) is established which, up to now, has not been verified yet for other general non-Markov estimators. A similar theorem is shown for the bootstrap, yielding resampling-based inference methods for statistical functionals. Formulas of the involved covariance functions are presented in detail. Particular applications include the conditional expected length of stay in a specific state, given occupation of another state in the past, as well as the construction of time-simultaneous confidence bands for the transition probabilities. The expected lengths of stay in the two-sample liver cirrhosis data-set by Andersen et al. (1993) are compared and confidence intervals for their difference are constructed. With borderline significance and in comparison to the placebo group, the treatment group has an elevated expected length of stay in the healthy state given an earlier disease state occupation. In contrast, the Aalen-Johansen estimator-based confidence interval, which relies on a Markov assumption, leads to a drastically different conclusion. Also, graphical illustrations of confidence bands for the transition probabilities demonstrate the biasedness of the Aalen-Johansen estimator in this data example. The reliability of these results is assessed in a simulation study.

Figures

Tables

Table 1. Table 1 : Empirical bias and coverage of nominal 95% confidence intervals for simulated datasets

	Bias		Coverage %
$n$	AJ	NM	Wald (AJ)	Wald (NM)	Naive bootstrap	Bootstrap-t
50	-0.231	0.181	92.09	91.01	94.23	90.85
100	-0.268	0.085	89.95	94.54	94.08	94.43
150	-0.279	0.063	88.14	95.24	94.21	95.00
200	-0.297	0.033	86.34	95.56	94.54	95.43

Table 2. Table 2 : Empirical coverage of nominal 95% simultaneous confidence bands for p 21 ​ ( 5 , t ) subscript 𝑝 21 5 𝑡 p_{21}(5,t) for t ∈ ( 6 , 7 ] 𝑡 6 7 t\in(6,7] . H-W = Hall-Wellner, EP = Equal precision

	Non-Markov			Aalen-Johansen
$n$	H-W	EP	Naive	EP
50	0.9996	0.9972	0.9226	0.9154
100	0.9896	0.9846	0.9238	0.9188
150	0.9604	0.9542	0.9298	0.9230
200	0.9496	0.9496	0.9420	0.9124

Table 3. Table 3 : Estimated difference in conditional length-of-stay for liver cirrhosis data

Method	$\Delta e_{10}$	95% CI
Wald	375.3	(-2.1, 752.8)
Naive bootstrap	375.3	(14.7, 740.0)
Bootstrap-t	375.3	(-5.1, 747.3)
AJ (Wald)	11.9	(-217.3, 241.1)

Equations

$$P_{\mathcal{I}\mathcal{J}}(s,t)=P(X(t)\in \mathcal{J}\mid X(s)\in \mathcal{I})$$

$$P_{\mathcal{I}\mathcal{J}}(s,t)=$$

$$+P(Y(t)=1\mid Y(s)=1,X(s)\in \mathcal{I})$$

$$\times P(X(t)\in \mathcal{J}\mid Y(t)=1,X(s)\in \mathcal{I})$$

$$=$$

$${\hat{p}}_{\mathcal{J}\mid \mathcal{I}}(t)=\frac{{\hat{p}}_{\mathcal{J}\mathcal{I}}(t)}{{\hat{p}}_{\mathcal{I}}(t)}=\frac{n^{-1}{\sum }_{\ell =1}^{n}1\{X_{\ell }(s)\in \mathcal{I},X_{\ell }(t)\in \mathcal{J},Y_{\ell }(t)=1\}}{n^{-1}{\sum }_{\ell =1}^{n}1\{X_{\ell }(s)\in \mathcal{I},Y_{\ell }(t)=1\}},$$

$${\hat{P}}_{\mathcal{I}\mathcal{J}}(s,t)={\hat{F}}_1(t)+{\hat{F}}_0(t){\hat{p}}_{\mathcal{J}|\mathcal{I}}(t).$$

$$\displaystyle\lim_{\delta\downarrow 0}\sup_{t-t_{1}=\delta}\Big{|}$$

$$\displaystyle\lim_{\delta\downarrow 0}\sup_{\begin{subarray}{c}t_{2}-t=\delta\\ t_{1}<t\end{subarray}}\Big{|}$$

$$\sqrt{n}({\hat{F}}_0-F_0,{\hat{F}}_1-F_1,{\hat{p}}_{\mathcal{J}|\mathcal{I}}-p_{\mathcal{J}|\mathcal{I}})\to (L_0,L_1,G_{\mathcal{J}|\mathcal{I}})$$

$$\left(\begin{array}{ccc}{\Sigma }_{00}& {\Sigma }_{01}& {\Sigma }_{0G}\\ {\Sigma }_{10}& {\Sigma }_{11}& {\Sigma }_{1G}\\ {\Sigma }_{G0}& {\Sigma }_{G1}& {\Sigma }_{GG}\end{array}\right):[s,\tau {]}^2\to {\mathbb{R}}^{3\times 3}.$$

$$\sqrt{n}({\hat{P}}_{\mathcal{I}\mathcal{J}}(s,\cdot )-P_{\mathcal{I}\mathcal{J}}(s,\cdot ))\to U$$

$$(r,t)$$

$$+p_{\mathcal{J}|\mathcal{I}}(r){\Sigma }_{01}(r,t)+p_{\mathcal{J}|\mathcal{I}}(t){\Sigma }_{10}(r,t)$$

$$+p_{\mathcal{J}|\mathcal{I}}(r\vee t)F_0(r\wedge t){\Sigma }_{0G}(r\vee t,r\wedge t)+F_0(r\wedge t){\Sigma }_{1G}(r\vee t,r\wedge t).$$

$$P_{\mathcal{I}\mathcal{J}}^{\ast }(s,t)=F_1^{\ast }(t)+F_0^{\ast }(t)p_{\mathcal{J}|\mathcal{I}}^{\ast }(t),t\in [s,\tau ].$$

$$\sqrt{n}(P_{\mathcal{I}\mathcal{J}}^{\ast }(s,\cdot )-{\hat{P}}_{\mathcal{I}\mathcal{J}}(s,\cdot ))\to U$$

$$e_{\mathcal{I}\mathcal{J}}(s,\tau )$$

$$={\int }_s^{\tau }P(X(u)\in \mathcal{J}|X(s)\in \mathcal{I})\mathrm{d}u;$$

$$H_{\mathcal{J}\mid \mathcal{I}}:e_{\mathcal{I}\mathcal{J}}^{(1)}(s,\tau )\ge e_{\mathcal{I}\mathcal{J}}^{(2)}(s,\tau )\text{versus}{K}_{\mathcal{J}\mid \mathcal{I}}:e_{\mathcal{I}\mathcal{J}}^{(1)}(s,\tau )<e_{\mathcal{I}\mathcal{J}}^{(2)}(s,\tau ).$$

$$\sqrt{\frac{n_1{n}_2}{n_1+n_2}}\frac{{\hat{e}}_{\mathcal{I}\mathcal{J}}^{(1)}(s,\tau )-{\hat{e}}_{\mathcal{I}\mathcal{J}}^{(2)}(s,\tau )}{{\hat{\sigma }}_{\mathcal{I}\mathcal{J}}}\to N(0,1).$$

$${\sigma }_{\mathcal{I}\mathcal{J}}^2=(1-\lambda ){\int }_s^{\tau }{\int }_s^{\tau }{\Gamma }_{\mathcal{I}\mathcal{J}}^{(1)}(u,v)\mathrm{d}u\mathrm{d}v+\lambda {\int }_s^{\tau }{\int }_s^{\tau }{\Gamma }_{\mathcal{I}\mathcal{J}}^{(2)}(u,v)\mathrm{d}u\mathrm{d}v$$

$$B_n:t⟼\sqrt{n}\cdot w(t)\cdot (\varphi ({\hat{P}}_{\mathcal{I}\mathcal{J}}(s,t))-\varphi (P_{\mathcal{I}\mathcal{J}}(s,t))),t\in [t_1,t_2],$$

$$w_{21}(t)$$

$$\displaystyle=-\log(\widehat{P}_{\mathcal{I}\mathcal{J}}(s,t))/\Big{[}1+\frac{\widehat{\sigma}^{2}_{\mathcal{I}\mathcal{J}}(t)}{\widehat{P}_{\mathcal{I}\mathcal{J}}^{2}(s,t)}\Big{]}^{-1},$$

$$\text{or}{w}_{22}(t)$$

$$\displaystyle=-\log(1-\widehat{P}_{\mathcal{I}\mathcal{J}}(s,t))/\Big{[}1+\frac{\widehat{\sigma}^{2}_{\mathcal{I}\mathcal{J}}(t)}{(1-\widehat{P}_{\mathcal{I}\mathcal{J}}(s,t))^{2}}\Big{]}^{-1}.$$

$$C{B}_{1k}$$

$$=[{\hat{P}}_{\mathcal{I}\mathcal{J}}(s,t{)}^{\exp (\mp q_{1-\alpha }^{\ast }/(\sqrt{n}{w}_{1k}(t)))}{]}_{t_1\le t\le t_2}$$

$$\text{and}C{B}_{2k}$$

$$=[1-(1-{\hat{P}}_{\mathcal{I}\mathcal{J}}(s,t){)}^{\exp (\mp q_{1-\alpha }^{\ast }/(\sqrt{n}{w}_{2k}(t)))}{]}_{t_1\le t\le t_2}$$

$${\lambda }_{01}(t)=\{\begin{array}{ll}0.3& \text{mbox}ifX(4)=1,t>4\\ 0.6& \text{mbox}otherwise.\end{array}$$

$$\displaystyle\int_{s}^{\tau}\int_{s}^{\tau}\Big{(}$$

$$\displaystyle+p_{0\mid 1}(u\vee v)F_{0}(u\wedge v)\Sigma_{0G}(u\wedge v,u\vee v)\Big{)}\mathrm{d}u\mathrm{d}v.$$

$$\sum _{j=0}^2\sum _{k\ne j}{\int }_s^{r\wedge t}{P}_{1j}^2(s,u)(P_{k0}(u,t)-P_{j0}(u,t))(P_{k0}(u,r)-P_{j0}(u,r))\frac{{\alpha }_{jk}(u)}{y_j(u)}\mathrm{d}u;$$

$${\Sigma }_{\mathcal{J}\mathcal{I}}(u,v)$$

Full text

In this article, weak convergence of the general non-Markov state transition probability estimator by Titman (2015) is established which, up to now, has not been verified yet for other general non-Markov estimators. A similar theorem is shown for the bootstrap, yielding resampling-based inference methods for statistical functionals. Formulas of the involved covariance functions are presented in detail. Particular applications include the conditional expected length of stay in a specific state, given occupation of another state in the past, as well as the construction of time-simultaneous confidence bands for the transition probabilities. The expected lengths of stay in the two-sample liver cirrhosis data-set by Andersen et al. (1993) are compared and confidence intervals for their difference are constructed. With borderline significance and in comparison to the placebo group, the treatment group has an elevated expected length of stay in the healthy state given an earlier disease state occupation. In contrast, the Aalen–Johansen estimator-based confidence interval, which relies on a Markov assumption, leads to a drastically different conclusion. Also, graphical illustrations of confidence bands for the transition probabilities demonstrate the biasedness of the Aalen-Johansen estimator in this data example. The reliability of these results is assessed in a simulation study.

Keywords: Aalen–Johansen estimator; conditional expected length of stay; confidence bands; Markov assumption; multi-state model; right-censoring; weak convergence.

1 Introduction

Transition probabilities are essential quantities in survival analytic examinations of multi-state models. Under independent right-censoring, for instance, the Aalen–Johansen estimator assesses these probabilities optimally in nonparametric Markovian multi-state models; cf. Aalen (1978) and Section IV.4 in Andersen et al. (1993) for more details. The Markov assumption is crucial for the validity of the Aalen–Johansen estimator. However, in real world applications, it may be unrealistic to postulate such a structure: given the present patient’s state, future developments are often not independent of a severe past illness history. There have been several attempts to circumvent the Markov assumption: State occupation probabilities in general models have been estimated by Datta and Satten (2001), Datta and Satten (2002), and Glidden (2002) using the Aalen–Johansen estimator. Pepe et al. (1991) and Pepe (1991) used a combination of two Kaplan–Meier estimators to assess the relevant transition probability in an illness-death model with recovery. In illness-death models without recovery, Meira-Machado et al. (2006) used Kaplan–Meier-based techniques and consecutive latent transition times to estimate transition probabilities. More efficient variants of these estimators have been developed in de Uña-Álvarez and Meira-Machado (2015). A Kendall’s $\tau$-based test of the Markov assumption in this progressive illness-death model was derived in Rodríguez-Girondo and Uña-Álvarez (2012). A competing risks-based estimator was proposed by Allignol et al. (2014) which relies on strict censoring assumptions. Eventually, Titman (2015) and Putter and Spitoni (2018) developed transition probability estimators in general non-Markov models by using different transition probability decompositions and utilizations of the Aalen–Johansen estimator. However, weak convergence properties of these estimators as elements of càdlàg function spaces have not been analyzed yet.

In the present paper, we focus on the non-Markov state transition probability estimator by Titman (2015) and its weak convergence for the following reasons: the estimator has a simple, but intuitive structure and it allows estimation of general transition probabilities between sets of states rather than single states. Finally, the bootstrap is shown to correctly reproduce the weak limit distribution. Throughout the article, let $(\Omega,\mathcal{A},P)$ be the underlying probability space. To introduce the Titman (2015) estimator, we consider an independently right-censored multi-state model with $R+1\in\mathbb{N}$ different states. That is, the data consists of independent copies $X_{1},\dots,X_{n}$ of a multi-state process which are indexed by time $t\geq s$ and which occupy the states $0,1,\dots,R$. The right-censoring is modeled by independent random variables $C_{1},\dots,C_{n}$ after which further observation of the corresponding multi-state processes is not possible any more. Let $\mathcal{I}$ and $\mathcal{J}$ be subsets of $\{0,1,\dots,R\}$ for which the transition probabilities

[TABLE]

shall be estimated. Here and below, irrelevant subscripts are removed for notational convenience. Introduce the set of states $\mathcal{A}_{\mathcal{J}}\subset\mathcal{J}$, that implies sure occupancy of $\mathcal{J}$ at all later points in time, and the set of states $\mathcal{R}_{\mathcal{J}}\subset\{0,1,\dots,R\}\setminus\mathcal{J}$ which prevents occupancy of $\mathcal{J}$ at all later points of time. Individual competing risks processes, $Z_{i},i=1,\dots,n$, are used in an intermediate estimation step. $Z_{i}(t)=0$ indicates whether the multi-state process $i$ is at time $t$ in none of the absorbing subsets $\mathcal{A}_{\mathcal{J}}$ or $\mathcal{R}_{\mathcal{J}}$. On the other hand, $Z_{i}(t)=1$ holds if the transition to $\mathcal{A}_{\mathcal{J}}$ has been observed until time $t$ and similarly $Z_{i}(t)=2$ for the transition to $\mathcal{R}_{J}$. Let $Y_{i}(t),i=1,\dots,n$ denote the at-risk indicators of the competing risks processes and let $\delta_{i},i=1,\dots,n$ be the corresponding censoring indicators, i.e., $\delta_{i}=1$ if the transition of $Z_{i}$ is eventually observed and $\delta_{i}=0$ if a censoring comes first. The final data-set is $\chi=\{X_{i},Z_{i},Y_{i},\delta_{i}:i=1,\dots,n\}$ whereof, in case of $\delta_{i}=0$, the processes $X_{i}$ and $Z_{i}$ are observable only until censoring occurs. Hence, incorporating random right-censoring, the transition probability is decomposed into

[TABLE]

The first two components in this decomposition are estimated by the standard Kaplan–Meier estimator $\widehat{F}_{0}$ for the survival function and the Aalen–Johansen estimator $\widehat{F}_{1}$ for the cumulative incidence function of the sub-group of those individuals which were observed to occupy state $\mathcal{I}$ at time $s$. The remaining conditional probability is estimated using the empirical proportion

[TABLE]

where $1\{\cdot\}$ denotes the indicator function. Therefore, the overall transition probability estimator by Titman (2015) is given as

[TABLE]

This article shall establish the weak convergence of $\sqrt{n}(\widehat{P}_{\mathcal{I}\mathcal{J}}(s,\cdot)-P_{\mathcal{I}\mathcal{J}}(s,\cdot))$ and of a bootstrap variant thereof as $n\rightarrow\infty$. Finally, implications to confidence interval construction for the conditional expected length of stay and to time-simultaneous confidence bands for the transition probabilities are demonstrated. The proofs of all Theorems and Lemmata as well as detailed derivations of all asymptotic covariance functions are given in the appendix.

2 Main Results

We focus on the Titman (2015) estimator during $[s,\tau]$ where $\tau\in(s,\infty)$ is a terminal time satisfying $F_{0}(\tau)>0$ and $P(C>\tau)>0$. The large sample properties of $\widehat{P}_{\mathcal{I}\mathcal{J}}(s,\cdot)$ are derived from central limit theorems of its individual components. To obtain such for the fraction process $\widehat{p}_{\mathcal{J}|\mathcal{I}}$, we throughout assume the existence of bounded one- and two-step hazard functions for transitions into $\mathcal{J}$ and $\bar{\mathcal{J}}=\{0,1,\dots,R\}\setminus\mathcal{J}$, respectively, which are approached uniformly:

[TABLE]

Here, $t_{1},t,t_{2}\in(s,\tau]$. In contrast to the standard choice in the literature, we let the at-risk indicators $Y_{i},i=1,\dots,n$, that are involved in the estimator $\widehat{p}_{\mathcal{J}|\mathcal{I}}$, be right-continuous. The typical reason for using the left-continuous version is the applicability of martingale methods which do not matter in our analysis of the estimators. On top of that, this choice does not matter in the present article as all transition times are continuously distributed. The definitions of the Kaplan–Meier and the Aalen–Johansen estimators remain unchanged. We equip the càdlàg function space $D[s,\tau]$ with the supremum norm and products thereof with the maximum-supremum norm; see Chapter 3.9 in van der Vaart and Wellner (1996) for applications in survival analysis with this choice. In order to state our main theorem, it is convenient to provide an auxiliary lemma:

Lemma 1

Under conditions (1) to (4) and as $n\rightarrow\infty$, the following joint convergence holds in distribution on $D^{3}[s,\tau]$:

[TABLE]

Here, $L_{0}$, $L_{1}$, and $G_{\mathcal{J}|\mathcal{I}}$ are zero-mean Gaussian processes. Denote their variance-covariance functions, arranged in a convenient $(3\times 3)$-matrix, by

[TABLE]

Detailed formulas for these covariance functions are presented in the appendix. From the presentation therein it is apparent that the covariances can be estimated straightforwardly.

Define $r\wedge t=\min(r,t)$ and $r\vee t=\max(r,t)$. An application of the functional delta-method immediately yields our main theorem:

Theorem 1

Under conditions (1) to (4) and as $n\rightarrow\infty$, the following convergence holds in distribution on $D[s,\tau]$:

[TABLE]

where $U$ is a zero-mean Gaussian process with covariance function $\Gamma_{\mathcal{I}\mathcal{J}}$ given by

[TABLE]

To be able to use this result for time-simultaneous inference on $P_{\mathcal{I}\mathcal{J}}(s,\cdot)$ or simply for improving the accurateness of $P_{\mathcal{I}\mathcal{J}}(s,\cdot)$-based inference methods with a pivotal limit distribution, the classical bootstrap seems to be a useful choice of a resampling method. To this end, a bootstrap sample $\chi^{*}$ is obtained by randomly drawing quadruples $(X_{i},Z_{i},Y_{i},\delta_{i})$ $n$ times with replacement from the original sample $\chi$. This bootstrap sample is used to recalculate the transition probability estimator, resulting in

[TABLE]

For this bootstrap counterpart, a central limit theorem similar to the main Theorem 1 holds:

Theorem 2

Under conditions (1) to (4), as $n\rightarrow\infty$, and given $\chi$, the following conditional convergence in distribution on $D[s,\tau]$ holds in probability:

[TABLE]

where $U$ is the same Gaussian process as in Theorem 1.

That is, the bootstrap succeeds in correctly recovering the limit distribution of the transition probability estimate $\widehat{P}_{\mathcal{I}\mathcal{J}}(s,\cdot)$, enabling an approximation of its finite sample distribution.

3 Inference Procedures

3.1 Conditional Expected Length of Stay

A practical application of both main theorems is a two-sample comparison of expected lengths of stay in the state set of interest, $\mathcal{J}$, given occupation of a fixed state set $\mathcal{I}$ at time $s\in[0,\tau]$. Considering temporarily the above one-sample set-up, the conditional expected length of stay is

[TABLE]

see e.g. Grand and Putter (2016) for a pseudo-observation regression treatment.

An R package for deriving the change in lengths of stay, that relies on the Aalen–Johansen estimator, is described in Wangler et al. (2006). For further literature on lengths of stay within multi-state models, see the articles cited therein. In the special case of a simple survival model, the only expected length of stay of interest is the mean residual lifetime; see e.g. Meilijson (1972). Here, however, the Markov assumption shall not be imposed such that estimation of the length of stay is based on the Titman (2015) estimator.

Sample-specific quantities are furnished with a superscript ${}^{(\ell)},\ell=1,2$. For instance, $\widehat{P}^{(1)}_{{\mathcal{I}\mathcal{J}}}(s,t)$ is the first sample’s transition probability estimate to switch into some state of $\mathcal{J}$ at time $t$ given that a state of $\mathcal{I}$ is occupied at time $s$. The sample sizes $n_{1}$ and $n_{2}$ may be different. We would like to test the null hypothesis

[TABLE]

If $\mathcal{J}$ is a favourable set of states, then the rejection of $H_{\mathcal{J}\mid\mathcal{I}}$ is an indication of the superiority of the second treatment over the first. A reasonable estimator of ${e}^{(\ell)}_{\mathcal{I}\mathcal{J}}(s,\tau)$ is $\widehat{e}_{\mathcal{I}\mathcal{J}}^{(\ell)}(s,\tau)=\int_{s}^{\tau}\widehat{P}_{\mathcal{I}\mathcal{J}}^{(\ell)}(s,u)\mathrm{d}u,$  $\ell=1,2,$. Note that, even under the boundary $\partial H_{\mathcal{J}|\mathcal{I}}:{e}^{(1)}_{{\mathcal{I}\mathcal{J}}}(s,\tau)={e}^{(2)}_{{\mathcal{I}\mathcal{J}}}(s,\tau)$, the underlying transition probabilities $P^{(1)}_{\mathcal{I}\mathcal{J}}$ and $P^{(2)}_{\mathcal{I}\mathcal{J}}$ are allowed to differ. Hence, we do not rely on the restrictive null hypothesis $\tilde{H}:P^{(1)}_{\mathcal{I}\mathcal{J}}\equiv P^{(2)}_{\mathcal{I}\mathcal{J}}$ but rather on the hypothesis $H_{\mathcal{J}\mid\mathcal{I}}$ of actual interest.

Theorem 3

Under conditions (1) to (4) and if  $0<\liminf n_{1}/n_{2}\leq\limsup n_{1}/n_{2}$ $<\infty$ as $n\rightarrow\infty$, the following convergence in distribution holds:

[TABLE]

Here, $\widehat{\sigma}^{2}_{\mathcal{I}\mathcal{J}}$ is any consistent estimator, e.g., found via plug-in, of the asymptotic variance

[TABLE]

along subsequences of $n_{1}/(n_{1}+n_{2})\rightarrow\lambda\in(0,1)$.

Replacing all estimators in the previous theorem by their bootstrap counterparts, where both underlying bootstrap samples $\chi^{(1)*}$ and $\chi^{(2)*}$ are obtained by $n_{1}$ and $n_{2}$ times independently and separately drawing with replacement from $\chi^{(1)}$ and $\chi^{(2)}$, respectively, a similar convergence in conditional distribution given $\chi^{(1)}\cup\chi^{(2)}$ holds in probability. This result may be used for bootstrap-based and asymptotically exact confidence intervals for the difference in the conditional expected lengths of stay, ${e}_{{\mathcal{I}\mathcal{J}}}^{(1)}(s,\tau)-{e}_{{\mathcal{I}\mathcal{J}}}^{(2)}(s,\tau)$; see Sections 4 and 5 below.

3.2 Simultaneous Confidence Bands

Even though the bootstrap procedure is useful but not strictly necessary for inference on the expected length of stay, it plays an essential role in time-simultaneous inference on the transition probabilities $t\mapsto P_{\mathcal{I}\mathcal{J}}(s,t)$. This is due to the unknown stochastic behaviour of the limit process $U$ which, again in the context of creating confidence bands, is also the reason for the inevitableness of resampling procedures for Aalen-Johansen estimators even if the Markov assumption is fulfilled; see Bluhmki et al. (2018) for theoretical justifications and their practical performance.

For the derivation of reasonable time-simultaneous $1-\alpha\in(0,1)$ confidence bands for $P_{\mathcal{I}\mathcal{J}}(s,\cdot)$ on a time interval $[t_{1},t_{2}]\subset[s,\tau]$ based on $\widehat{P}_{\mathcal{I}\mathcal{J}}(s,\cdot)$ in combination with the bootstrap, we focus on the process

[TABLE]

where a transformation $\phi$ will ensure bands within the probability bounds $0\%$ and $100\%$, and a suitable multiplicative weight function $w$ will stabilize the bands for small samples. In particular, we use one of the log-log transformations $\phi_{1}(p)=\log(-\log(p))$ or $\phi_{2}(p)=\log(-\log(1-p))$ (Lin, 1997), depending on whether we expect $P_{\mathcal{I}\mathcal{J}}(s,t_{1})$ to be closer to 1 or 0. As weight function, for transformation $k\in\{1,2\}$, we choose $w_{1k}(t)=[\phi_{k}^{\prime}(\widehat{P}_{\mathcal{I}\mathcal{J}}(s,t))\cdot\widehat{\sigma}_{\mathcal{I}\mathcal{J}}(t)]^{-1}$,

[TABLE]

With these choices, the resulting confidence bands will correspond to the classical equal precision (EP) bands for $w=w_{1k}$ or the Hall-Wellner (H-W) bands for $w=w_{2k}$ if $\mathcal{J}$ is an absorbing subset of states; cf. p. 265ff. in Andersen et al. (1993) for the survival case and $k=1$ or Lin (1997) in the presence of competing risks and $k=2$.

The confidence bands are found by solving for $P_{\mathcal{I}\mathcal{J}}(s,t)$ in the probability $P(\sup_{t\in[t_{1},t_{2}]}|B_{n}(t)|\leq q_{1-\alpha})=1-\alpha$. Here, the value of $q_{1-\alpha}$ is approximated by the $(1-\alpha)$-quantile $q_{1-\alpha}^{*}$ of the conditional distribution of the supremum of the bootstrap version $B_{n}^{*}(t)=\sqrt{n}\cdot w^{*}(t)\cdot(\phi(P_{\mathcal{I}\mathcal{J}}^{*}(s,t))-\phi(\widehat{P}_{\mathcal{I}\mathcal{J}}(s,t)))$ of $B_{n}$ given the data, where in the definition of $w^{*}$ all estimators $\widehat{P}_{\mathcal{I}\mathcal{J}}(s,t)$ and $\widehat{\sigma}_{\mathcal{I}\mathcal{J}}(t)$ are replaced with their bootstrap counterparts $P_{\mathcal{I}\mathcal{J}}^{*}(s,t)$ and $\sigma^{*}_{\mathcal{I}\mathcal{J}}(t)$.

Finally, the resulting $1-\alpha$ confidence bands for $t\mapsto P_{\mathcal{I}\mathcal{J}}(s,t)$ are as follows:

[TABLE]

The performance of these confidence bands as well as the confidence band obtained from not transforming or weighting the transition probability estimator are assessed in a simulation study in the subsequent sections.

4 Simulation Study

4.1 Conditional Expected Length of Stay

To assess the performance of the conditional length of stay estimators for finite sample sizes a small simulation study is conducted. Data are simulated from a pathological non-Markov model in which the subsequent dynamics of the process depend on the state occupied at $t=4$. A three-state illness-death model with recovery is used where the state [math] means “healthy”, $1$ “ill”, and $2$ “dead”. The transition intensities are $\lambda_{02}(t)=0.02,\lambda_{10}(t)=0.3,\lambda_{12}(t)=0.1$, while

[TABLE]

In the last part of the appendix we prove that the conditions (1)–(4) are satisfied in the present set-up. Subjects are independently right-censored via an exponential distribution with rate 0.04. We focus on the expected length of stay in the health state conditional on an earlier illness at time $s=5$, i.e. $\mathcal{I}=\{1\}$ and $\mathcal{J}=\{0\}$. For each simulated dataset, $e_{10}(5,30)$ is estimated and 95% confidence intervals are constructed via three methods; a Wald interval using the plug-in variance estimator $\widehat{\sigma}_{10}$ of ${\sigma}_{10}$, in which the canonical estimators of $p_{0\mid 1}$, $F_{0}$, $\Sigma_{00}$, $\Sigma_{0G}$, and $\Sigma_{GG}$ are chosen, a classical bootstrap and a bootstrap-t procedure using the bootstrap-version of $\widehat{\sigma}_{10}$ to studentize the bootstrap samples. For each bootstrap, $B=1000$ samples are generated.

In the present three-state model, the Titman (2015) transition probability estimator reduces to {\color[rgb]{0,0,0}\widehat{P}_{10}(5,t)=}\widehat{P}_{\mathcal{I}\mathcal{J}}(5,t)=\widehat{F}_{0}(t)\widehat{p}_{\mathcal{J}|\mathcal{I}}(t){\color[rgb]{0,0,0}=\widehat{F}_{0}(t)\widehat{p}_{0\mid 1}(t)} because the healthy state is non-absorbing. Thus, the asymptotic variance of each sample-specific conditional expected length of stay estimator is

[TABLE]

As a competing method, Aalen–Johansen estimator-based Wald-type confidence intervals are constructed for the same expected length of stay. Standard arguments yield that the asymptotic covariance function of the normalized Aalen–Johansen estimator for the $1\rightarrow 0$ transition is $(r,t)\mapsto$

[TABLE]

cf. Section IV.4.2 in Andersen et al. (1993) for the corresponding particular variance formula. The expected length of stay asymptotic variance again results from a double integral of this covariance function with integration range $r,t\in(s,\tau]$.

We consider several different scenarios for the total sample size, $n=50,100,$ $150,$ and $200$. Note that the non-Markov estimator only uses the subgroup of subjects satisfying $C_{i}>5,X_{i}(t)=1$ which corresponds to only $0.367n$ subjects, on average. For each of the scenarios, 10000 datasets are generated. As a result the empirical coverage percentages have approximate Monte-Carlo standard error of around 0.2.

The simulation results are shown in Table 1. Both the Aalen-Johansen and non-Markov estimator have substantial bias when $n=50$ leading to under coverage of all the nominal 95% confidence intervals. For larger $n$, both of the bootstrap confidence intervals and the Wald-based interval for the non-Markov estimator give close to nominal coverage. The coverage of the Aalen-Johansen based Wald confidence interval deteriorates with increasing $n$ due to the estimator being inherently biased.

4.2 Simultaneous Confidence Bands

To investigate the coverage probabilities of simultaneous confidence bands introduced in Section 3.2, datasets are simulated under the same non-Markov model as above, again considering total sample sizes of $n=50,100,150,$ and $200$ patients. 95% simultaneous confidence bands are constructed for $p_{21}(5,t)$ for $t\in(6,7].$ For the non-Markov estimator, in addition to the Hall-Wellner and equal precision bands using the transformation $\phi_{2}(p)$, a naive confidence band based on a constant weight function, $w(t)\equiv 1$, and an identity transformation, $\phi(p)=p$, is also constructed. In addition, for the Aalen-Johansen estimates, EP bands based on $\phi_{2}(p)$, are constructed via the wild bootstrap using the R code which has been made available in the supplement to Bluhmki et al. (2018). In all cases, for each bootstrap $B=1000$ samples are generated.

Table 2 gives the empirical coverage probabilities of the simultaneous confidence bands based on 5000 simulated datasets. For the bands based on the Non-Markov estimator, the H-W and EP bands tend to over-cover quite markedly for small sample sizes, while the naive bands under-cover. Adequate coverage is achieved for the H-W and EP bands by $n=200$. Since the Aalen-Johansen estimator is biased for this scenario we would expect the coverage of the EP Aalen-Johansen bands to deteriorate with increasing sample size. However, it appears this is counteracted by the estimated quantiles, $q^{*}_{1-\alpha}$ increasing with $n$. As a consequence, the coverage initially grows with increasing $n$ before beginning to decrease.

5 Application to the Liver Cirrhosis Data-Set

As an illustrative example we consider the liver cirrhosis data set introduced by Andersen et al. (1993) and also analyzed by Titman (2015). Patients were randomized to either a treatment of prednisone (251 patients) or a placebo (237 patients). A three-state illness-death model with recovery is assumed where the healthy and illness states correspond to normal and elevated levels of prothrombin, respectively. A potential measure of the effectiveness of prednisone as a treatment is the expected length-of-stay in the normal prothrombin level state, from a defined starting point. Specifically we consider $\Delta e_{10}=e_{10}^{(1)}(s,\tau)-e_{10}^{(2)}(s,\tau)$, the difference in conditional expected length-of-stay in normal prothrombin levels, given the patient is alive and with abnormal prothrombin at time $s$, for the prednisone and placebo groups, corresponding to $\ell=1$ and $\ell=2$, respectively. Titman (2015) observed an apparent treatment difference with respect to the transition probabilities with respect to {\color[rgb]{0,0,0}s}=1000 days post-randomization. We consider the expected length-of-stay in state 0 up to time $\tau=3000$ days.

In this case, $n_{1}=26$ patients in the prednisone group and $n_{2}=35$ patients in the placebo group meet the condition that $X_{i}(1000)=1$ and $Y_{i}(1000)=1$. Bootstrap based confidence intervals are constructed using $1000$ bootstrap samples within each group.

Table 3 shows the three constructed 95% confidence intervals, which are broadly similar, although the naive bootstrap interval excludes 0, whereas the Wald and bootstrap-t intervals do not. The Aalen-Johansen based estimate and a Wald confidence interval is also given. It is seen that there is a dramatic difference between the non-Markov and Aalen-Johansen estimates, with the latter indicating no treatment difference. Potentially, the apparent non-Markov behaviour in the data may be due to patient heterogeneity within the treatment groups.

To illustrate the construction of simultaneous confidence bands we construct 95% confidence bands for $P_{21}(500,t)$, in the interval $t\in(750,1250]$. Since it was seen in Section 4 that reasonable sample size is required for good coverage an earlier start time is used here than for the estimate of the expected length of stay to ensure sufficient numbers of patients are under observation. Since we expect $P_{21}(500,t)<0.5$, we chose the confidence bands using the transformation $\phi_{2}(p)=\log(-\log(1-p)).$

Figure 1 shows the confidence bands for the placebo and prednisone groups using Hall-Wellner and equal precision intervals using the non-Markov estimator and using the equal precision intervals for the Aalen-Johansen estimator.

6 Discussion

To the best of our knowledge, our paper is the first to prove weak convergence properties of a general transition probability estimator in independently right-censored non-Markov multi-state models while also providing explicit formulas of the asymptotic covariance functions. Similar proofs have been given for the classical bootstrap, which finally allows the utilization of the Titman (2015) estimator in time-simultaneous inference procedures. Unfortunately, when focusing on inference on the expected length of stay for single time points, both bootstrap methods applied in Section 4 yielded a greater deviance of the nominal coverage level than the simple studentization method based on the asymptotic quantiles of the standard normal distribution. Even though the simulated small sample coverage probabilities were more satisfying, other resampling procedures might improve the performance in both, small and large samples. For example, the martingale residual multiplier methods of Lin (1997), Beyersmann et al. (2013), or Dobler et al. (2017) are the method of choice in incomplete competing risks data and they may be adapted for the present context in a future paper.

On the other hand, the time-simultaneous confidence bands derived in this paper cannot rely on a purely asymptotic quantile finding approach. Instead, they truly require a resampling method such as the bootstrap method developed here. Furthermore, the simulation results in Section 4.2 and the application to the liver cirrhosis data set in Section 5 demonstrate their practical usefulness: not only did they reveal reliable coverage probabilities for relatively small sample sizes already. They were also only slightly wider than the confidence bands based on the Aalen-Johansen estimator for which there is no guarentee of applicability.

Another field of possible future interest is the extension of the present methodology to quality-adjusted sums of expected lengths of stay in a specific sub-set of states; cf. Williams (1985) who proposed such quality-adjusted life years for operation considerations. For instance, the expected length of stay in a healthy state could be weighted with a factor $p\in(0.5,1)$ whereas the length of stay in the illness state obtains the weight $1-p$. The specific choice of $p$ may be obtained from additional information on, e.g., the pain scores of patients in the illness state and may be chosen individually by the patients themselves. This way, a comparison of treatments could be achieved that more realistically accounts for the circumstances of a specific disease. However, one should bear in mind the controversial debate on quality-adjusted life years as articulated by e.g. Harris (1987); see also the discussion of the paper by Cox et al. (1992). Therein, G. W. Torrance appreciates the usefulness of quality-adjusted life years for resource allocation if combined with other measures and advices:

He states that quality-adjusted life years “were never intended for clinical decision-making - they were developed for use in resource allocation.” They should also be reported together with other important measures for reaching a decision. Finally, one should add “a thorough sensitivity analysis and thoughtful discussion, including caveats” to obtain an informative quality-adjusted life years analysis.

Finally, reconsider the data analysis in Section 5 where the Aalen–Johansen estimator yielded a difference of expected lengths of stay which drastically deviates from the same quantity based on the Titman (2015) estimator. This might be considered as a strong evidence against the Markov assumption. If, on the other hand, a formal test of Markovianity does not find such evidence, the resulting estimated lengths of stay might be comparable as well. Even more can be gained by means of such tests: as apparent from the mentioned liver cirrhosis data example, the confidence intervals and simultaneous confidence bands based on the Aalen–Johansen estimator are much narrower, resulting from the efficiency of this estimator in comparison to non-Markov transition probability estimates. Equivalently, test procedures based on the classical estimator are more powerful in detecting deviances from null hypotheses. Therefore, formal tests for the applicability of the classical Aalen–Johansen estimator shall be developed in future articles.

Acknowledgements

The authors would like to thank Markus Pauly for helpful discussions. Dennis Dobler was supported by a DFG (German Research Foundation) grant.

Appendix A Asymptotic covariance functions and comments on their estimation

We begin by stating all covariance functions which are involved in the limit theorems. In particular, introducing the following covariance functions will prove to be useful:

[TABLE]

With these abbreviations, we conveniently introduce the asymptotic covariance functions which appear in Lemma 1 of the manuscript. Their derivation in detail will be done in the next section.

[TABLE]

The asymptotic variance-covariance functions of the Kaplan–Meier and the Aalen–Johansen estimator for different time points are well-known and therefore not listed. The derivation of $cov(L_{0}(r),G_{\mathcal{J}|\mathcal{I}}(u))$ and $cov(L_{1}(t),G_{\mathcal{J}|\mathcal{I}}(u))$ will be made explicit in Section F below. In the display above, we used $P_{\mathcal{I}\mathcal{J}k}(s,u,t)=P(Z(t)=k,X(u)\in\mathcal{J}\mid X(s)\in\mathcal{I}),k=1,2$. Note that these functions may be decomposed into

[TABLE]

of which the first probability is again estimable by Kaplan–Meier or Aalen–Johansen-type estimators based on those individuals that are observed to satisfy X(u)\in{\color[rgb]{0,0,0}\mathcal{J}} and $X(s)\in\mathcal{I}$. Hence, consistent plug-in estimates of the above covariance function are found easily. However, the number of individuals satisfying that X(u)\in{\color[rgb]{0,0,0}\mathcal{J}} and that $X(s)\in\mathcal{I}$ is observed might be quite small in practical applications. This could lead to a large variance of the estimators of $P_{\mathcal{I}\mathcal{J}k}(s,u,t)$.

Appendix B Proof of Lemma 1

The large sample properties concerning the conditional Kaplan–Meier estimator $\widehat{F}_{0}$ and the conditional Aalen–Johansen estimator $\widehat{F}_{1}$ follow from standard theory. Their variance-covariance functions are obtained in the same way. The asymptotic covariance functions between $\widehat{p}_{\mathcal{J}|\mathcal{I}}$ and $\widehat{F}_{0}$ as well as $\widehat{F}_{1}$ are derived in Appendix F below. For technical reasons, we consider weak convergence on the càdlàg space $D[s,\tau]$ equipped with the supremum norm instead of the usual Skorohod metric. In case of continuous limit sample paths, weak convergence on the Skorohod and the supremum-normed càdlàg space are equivalent; cf. Gill (1989), p. 110 in Andersen et al. (1993), or p. 137 in Pollard (1984). The advantage is that the functional delta-method is available on the supremum-normed space. Additionally, the general Aalen–Johansen estimator is known to have continuous limit distribution sample paths; cf. Theorem IV.4.2 in Andersen et al. (1993). However, this may not be true for the weak limit of $\sqrt{n}(\widehat{p}_{\mathcal{J}|\mathcal{I}}-p_{\mathcal{J}|\mathcal{I}})$ as $n\rightarrow\infty$ because the censoring distribution may have discrete components. This problem is solved by an application of a technique as used by Akritas (1986): Each discrete censoring components is distributed uniformly on adjacent, inserted small time intervals during which no state transition can occur. Instead, all future state transitions are shifted by the preceding interval lengths. Write $W_{\mathcal{J}\mathcal{I}}=\sqrt{n}(\widehat{p}_{\mathcal{J}\mathcal{I}}-p_{\mathcal{J}\mathcal{I}})$ and $W_{\mathcal{I}}=\sqrt{n}(\widehat{p}_{\mathcal{I}}-p_{\mathcal{I}})$, where $p_{{\mathcal{J}\mathcal{I}}}(t)=P(X(s)\in\mathcal{I},X(t)\in\mathcal{J},Y(t)=1)$ and $p_{\mathcal{I}}(t)=P(X(s)\in\mathcal{I},Y(t)=1)$. The resulting weak limit sample paths of the time transformed $\sqrt{n}(\widehat{p}_{\mathcal{I}}-p_{\mathcal{I}})$ and $\sqrt{n}(\widehat{p}_{\mathcal{J}\mathcal{I}}-p_{\mathcal{J}\mathcal{I}})$ will be continuous and the functional delta-method is applicable to obtain the corresponding result for $\widehat{p}_{\mathcal{J}|\mathcal{I}}$. Note that the values of none of these processes were modified outside the just inserted small time intervals. Finally, the projection functional, which projects a time-transformed process to the original time line, is continuous. Hence, the continuous mapping theorem eventually yields the desired central limit theorem for the process $\sqrt{n}(\widehat{p}_{\mathcal{J}|\mathcal{I}}-p_{\mathcal{J}|\mathcal{I}})$ as $n\rightarrow\infty$.

It remains to give a weak convergence proof of the finite-dimensional margins of the process $\sqrt{n}(\widehat{p}_{\mathcal{J}|\mathcal{I}}-p_{\mathcal{J}|\mathcal{I}})$ and to verify its tightness under the temporary assumption of continuously distributed censoring times. Combining both properties, weak convergence on $D[s,\tau]$ follows, as indicated above.

Finite-Dimensional Marginal Distributions

Weak convergence of the finite-dimensional margins is shown by applying the functional delta-method after utilizations of the central limit theorem. The weak convergences of all finite-dimensional marginal distributions are easily obtained by the multivariate central limit theorem: For $m\in\mathbb{N}$ points of time, $s\leq t_{1}<\dots<t_{m}\leq\tau$, consider $(W_{{\mathcal{J}\mathcal{I}}}(t_{1}),\dots,W_{{\mathcal{J}\mathcal{I}}}(t_{m}))^{\prime}\rightarrow(\tilde{G}_{{\mathcal{J}\mathcal{I}}}(t_{1}),\dots,\tilde{G}_{{\mathcal{J}\mathcal{I}}}(t_{m}))^{\prime}\sim N(0,\Sigma_{{\mathcal{J}\mathcal{I}}})$ in distribution, as $n\rightarrow\infty$, with $(j,k)$th covariance entry

[TABLE]

Similarly, $(W_{\mathcal{I}}(t_{1}),\dots,W_{\mathcal{I}}(t_{m}))^{\prime}\rightarrow(\tilde{G}_{\mathcal{I}}(t_{1}),\dots,\tilde{G}_{\mathcal{I}}(t_{m}))^{\prime}\sim N(0,\Sigma_{\mathcal{I}})$ in distribution, as $n\rightarrow\infty$, with $(j,k)$th covariance entry

[TABLE]

The convergence of both vectors even holds jointly, implying that the $(j,k)$th entry of the covariance matrix $\Omega_{\mathcal{J}\mathcal{I}}$ of both limit vectors is given by

[TABLE]

Denote by $D_{>0}[s,\tau]$ the subset of càdlàg functions with are positive and bounded away from zero, which is equipped with the supremum norm. Note that $\phi:D[0,\tau]\times D_{>0}[s,\tau]\rightarrow D[s,\tau],(f,g)\mapsto\frac{f}{g}$ defines a Hadamard-differentiable map with Hadamard-derivative $\phi^{\prime}_{(f,g)}(h_{1},h_{2})=\frac{h_{1}}{g}-h_{2}\frac{f}{g^{2}}$, with $h_{1},h_{2}\in D[s,\tau]$. This of course also implies Fréchet-differentiability on all finite-dimensional projections. Therefore, the multivariate delta-method yields the convergence in distribution

[TABLE]

as $n\rightarrow\infty$. The asymptotic covariance matrix $\Gamma_{{\mathcal{J}|\mathcal{I}}}$ has $(j,k)$th entry

[TABLE]

Tightness

We will show tightness on the Skorohod space $D[s,\tau]$ of càdlàg functions by showing that a modulus of continuity becomes small in probability; see Theorem 13.5 in Billingsley (1999). Remember that this property will imply weak convergence on the Skorohod space $D[s,\tau]$ which, due to the temporary continuity of the limit processes, is equivalent to weak convergence on the supremum-normed space.

Now, for the modulus of continuity, we consider for any small $\delta>0$ the expectation

[TABLE]

We start by rewriting $W_{\mathcal{I}}(t)-W_{\mathcal{I}}(t_{1})$ as

[TABLE]

Note that, in the product appearing in (8), the products of the above indicator functions with identical indices vanish for unequal time pairs $(t_{1},t)$ and $(t,t_{2})$. Additionally, due to the independence of all individuals, all terms with a single separate index $\ell$ vanish too. Hence, by an application of the Cauchy-Schwarz inequality, the expectation in (8) is bounded above by

[TABLE]

for all $n\geq 3$. Let $S_{C}(t)=\textnormal{pr}(C>t)$ denote the survival function of the censoring times and let $r_{1}(t)=P(Y(t)=1)\leq S_{C}(t)$ denote the (unconditional) probability to be still under risk at time $t$. By the temporarily continuous censoring distribution, $r_{1}$ is a continuous and non-decreasing function and thus it is suitable for an application of Theorem 13.5 in Billingsley (1999). Then the terms in the above display are not greater than

[TABLE]

Hence, tightness in the Skorohod space $D[s,\tau]$ of the process $\sqrt{n}(\widehat{p}_{\mathcal{I}}-p_{\mathcal{I}})$ affiliated with the denominator of $\widehat{p}_{{\mathcal{J}|\mathcal{I}}}$ follows.

It remains to prove tightness of the numerator process. At first, we rewrite for $t_{1}<t<t_{2}$

[TABLE]

For later calculations, we derive products of the $\zeta$-indicators for different points of time:

[TABLE]

Therefore, the expectation to be bounded for verifying tightness has the following structure:

[TABLE]

where the inequality follows from an application of the Cauchy-Schwarz inequality. The terms in the previous display are not greater than

[TABLE]

It remains to calculate and bound the single expectations.

[TABLE]

holds by the same reasons as above and a similar representation holds for the other second moments. By conditions (1) to (4) in the manuscript, there is a global constant $K>0$, that is independent of $t_{1},t,t_{2},n$, such that q_{\mathcal{I}01}(s,t_{1},t)\leq({\color[rgb]{0,0,0}\lambda_{\bar{\mathcal{J}}\mathcal{J}}(t)}+K)(t-t_{1}). Similarly, q_{\mathcal{I}10}(s,t_{1},t)\leq({\color[rgb]{0,0,0}\lambda_{\mathcal{J}\bar{\mathcal{J}}}(t)}+K)(t-t_{1}) as well as $q_{\mathcal{I}11}(s,t_{1},t)\leq r_{1}(t_{1})-r_{1}(t)$.

The first expectation (involving fourth moments) contains the following terms:

[TABLE]

For the expectation, first we obtain

[TABLE]

Similarly, we have $e_{2},e_{3},e_{31},e_{23},e_{32}\leq 0$. It remains to consider $e_{12},e_{13}$, and $e_{21}$.

[TABLE]

Similarly, $e_{12}\leq P(X_{1}(t_{1})\notin\mathcal{J},X_{1}(t)\in\mathcal{J},X_{1}(t_{2})\notin\mathcal{J})$ and $e_{13}\leq P(X_{1}(t_{1})\notin\mathcal{J},X_{1}(t)\in\mathcal{J},C_{1}\in(t,t_{2}])=P(X_{1}(t_{1})\notin\mathcal{J},X_{1}(t)\in\mathcal{J})P(C_{1}\in(t,t_{2}])$ due to independent right-censoring.

Assume that $e_{12},e_{13},e_{21}>0$. Otherwise, nothing needs to be shown. Thus,

[TABLE]

Due to condition (4), there is a constant $K>0$, which is independent of $t_{1},t,t_{2},n$, such that

[TABLE]

for a global constant $\tilde{K}>0$. Similarly $e_{12}\leq\tilde{K}^{2}(t_{2}-t_{1})^{2}$. Therefore, we conclude that

[TABLE]

for the non-decreasing and continuous function $H:x\mapsto x+1-S_{C}(x)$.

All in all, we conclude that

[TABLE]

where $F(x)=9\tilde{K}(x+1-r_{1}(x))$ defines a continuous and non-decreasing function, independently of $t_{1},t,t_{2},n$. Hence, tightness of the process $W_{{\mathcal{J}\mathcal{I}}}$ follows which in turn implies convergence in distribution on the Skorohod space $D[s,\tau]$. Finally, proceed as indicated in the beginning of this Appendix section.

Appendix C Proof of Theorem 1

The proof of Theorem consists of an application of the functional delta-method applied to the $D^{3}[s,\tau]\rightarrow D[s,\tau]$ functional $(f_{0},f_{1},p)\mapsto f_{1}+f_{0}p$ which is Hadamard-differentiable at $(F_{0},F_{1},p_{\mathcal{J}|\mathcal{I}})$ with derivative

[TABLE]

For the asymptotic covariance function, note that both terms $cov(L_{0}(t),G_{\mathcal{J}|\mathcal{I}}(r))$ and $cov(L_{1}(t),G_{\mathcal{J}|\mathcal{I}}(r))$ vanish for $r\geq t$.

Appendix D Proof of Theorem 2

First, the finite-dimensional, conditional marginal distributions of the estimator $\widehat{P}_{ij}^{*}$ given $\chi$ are shown to be consistent, as $n\rightarrow\infty$ for the corresponding marginal distributions of the limit process $L_{1}+L_{0}p_{{\mathcal{J}|\mathcal{I}}}+G_{{\mathcal{J}|\mathcal{I}}}F_{0}$. This is verified via empirical process theory: Let $s\leq t_{1}<\dots<t_{d}\leq\tau$ be finitely many time points. The finite-dimensional marginals of the estimator $\widehat{P}_{\mathcal{J}|\mathcal{I}}(s,\cdot)$ are Hadamard-differentiable functionals of the empirical process based on

[TABLE]

where $T_{\ell}=\inf\{t\geq s:Y_{\ell}(t)\neq 1\}\wedge\tau$. The empirical process is indexed by the family of functions

[TABLE]

where the particular functions are defined as

[TABLE]

Certainly, $\mathcal{F}$ defines a Donsker class. The derivation of the Kaplan–Meier estimator based on parts of this empirical process has been verified in Section 3.9 of van der Vaart and Wellner (1996). This similarly holds for the Aalen–Johansen estimator under competing risks by another, final application of the Wilcoxon functional and the use of the cause-specific Nelson–Aalen estimator; see Section F below for a more detailed derivation. The proportion-type estimator $\widehat{p}_{\mathcal{J}|\mathcal{I}}(t_{j}),j=1,\dots,d,$ is given by a simple Hadamard-differentiable functional of the empirical process indexed by the functions $f_{0},f_{j},g_{1t}$. Therefore, combining again all separate estimators, the functional delta-method is applicable to the bootstrap empirical process and the same Hadamard-differentiable functional as well (in outer probability); see Chapter 3.6 in van der Vaart and Wellner (1996) for details. It follows that the conditional finite-dimensional marginal distributions of $\sqrt{n}(P^{*}_{{\mathcal{I}\mathcal{J}}}(s,\cdot)-\widehat{P}_{{\mathcal{I}\mathcal{J}}}(s,\cdot))$ converge weakly to the limit distributions of the original normalized estimator.

Second, conditional tightness in probability needs to be verified. To this end, we utilize a variant of the tightness criterion in Theorem 13.5 in Billingsley (1999), in which the involved continuous, non-decreasing function $F$ may depend on the sample size $n$ and it need not be continuous, but only point-wise convergent and asymptotically continuous. See the Appendix A.1 in the PhD thesis by Dobler (2016) or the comment on p. 356 in Jacod and Shiryaev (2003) for details. As in the first part of this appendix, we assume without loss of generality that the censoring distribution is continuous. To verify tightness, we consider the conditional version of expectation that we bounded accordingly for verifying tightness of the original, normalized estimator in the first part of this appendix. All quantities with a superscript star are the obvious bootstrap counterparts of the quantities introduced earlier. We again focus on proving tightness for the bootstrap counterparts $W_{{\mathcal{J}\mathcal{I}}}^{*}=\sqrt{n}(p_{{\mathcal{J}\mathcal{I}}}^{*}-\widehat{p}_{{\mathcal{J}\mathcal{I}}})$ and $W_{\mathcal{I}}^{*}=\sqrt{n}(p_{\mathcal{I}}^{*}-\widehat{p}_{\mathcal{I}})$ of $W_{{\mathcal{J}\mathcal{I}}}=\sqrt{n}(\widehat{p}_{{\mathcal{J}\mathcal{I}}}-p_{{\mathcal{J}\mathcal{I}}})$ and $W_{\mathcal{I}}=\sqrt{n}(\widehat{p}_{\mathcal{I}}-p_{\mathcal{I}})$, respectively. Furthermore, the normalized, bootstrapped Kaplan–Meier estimator is well-known to be tight and the same property is easily seen for the normalized, bootstrapped Aalen–Johansen estimator under competing risks. Finally, we will apply the continuous mapping theorem in order to conclude asymptotic tightness of $\sqrt{n}(P_{{\mathcal{J}\mathcal{I}}}^{*}(s,\cdot)-\widehat{P}_{{\mathcal{J}\mathcal{I}}}(s,\cdot))$.

Let us abbreviate $W^{*}_{{\mathcal{J}\mathcal{I}}}(t)-W^{*}_{{\mathcal{J}\mathcal{I}}}(t_{1})=\frac{1}{\sqrt{n}}\sum_{\ell=1}^{n}(m_{\ell}-1)a_{\ell}$ and $W^{*}_{{\mathcal{J}\mathcal{I}}}(t_{2})-W^{*}_{{\mathcal{J}\mathcal{I}}}(t)=\frac{1}{\sqrt{n}}\sum_{\ell=1}^{n}(m_{\ell}-1)b_{\ell}$, where $(m_{1},\dots,m_{n})\sim Mult(n,\frac{1}{n},\dots,\frac{1}{n})$ is independent of $\chi$. Let $s\leq t_{1}\leq t\leq t_{2}\leq\tau$. Denote by $K>0$ a constant which may become larger along the subsequent inequalities, but which is independent of $n,t,t_{1},t_{2}$. Then, by the Cauchy-Schwarz inequality,

[TABLE]

As shown in e.g. Dobler and Pauly (2014), we have E\Big{(}\prod_{\ell=1}^{4}(m_{\ell}-1)\Big{)}=O(n^{-2}) and $E(m_{1}-1)^{2}(m_{2}-1)(m_{3}-1)=O(n^{-1})$. Furthermore, it is easy to show that $E(m_{1}-1)^{4},E(m_{1}-1)^{2}(m_{2}-1)^{2},E(m_{1}-1)^{3}(m_{2}-1)<\infty$. Therefore, the expectation in the previous display is bounded above by

[TABLE]

We apply the Cauchy-Schwarz inequality to $\sum_{\ell=1}^{n}a_{\ell}$ (and similarly for the $b_{k}$’s) which is bounded above by $(n\sum_{\ell=1}^{n}a_{\ell}^{2})^{1/2}$. Therefore, the terms in the previous display are bounded above by

[TABLE]

which we abbreviate by $K/n^{2}(AB+A^{3/2}B^{1/2}+A^{1/2}B^{3/2})$. Now, verifying the unconditional Billingsley criterion implies conditional tightness in probability: Suppose $\vartheta_{n}^{*}$ is a random element of $D[s,\tau]$ for which tightness in probability given some random element $\vartheta_{n}$ shall be shown and for which it hold that: For all $\eta,\varepsilon>0$ we have for the modulus of continuity $w^{\prime\prime}$, using the notation of Billingsley (1999), that

[TABLE]

Then, by the Markov inequality,

[TABLE]

Therefore, we need to find a suitable upper bound for the expectation of (9). By the Cauchy-Schwarz inequality we have the upper bound

[TABLE]

Here, we begin by looking for a bound of

[TABLE]

We note that each $a_{\ell}$ and $b_{\ell}$ are differences of indicator functions. In particular, we have that

[TABLE]

Therefore, we can bound

[TABLE]

A similar bound holds for $b_{\ell}^{2}$. Thus, we have for the product

[TABLE]

Similarly to the proof in Appendix B, we obtain the upper bound

[TABLE]

due to conditions (1) to (4) in the main article, where the function $H$ was defined in Appendix B.

Next, we find an upper bound for

[TABLE]

which holds for the same reasons as before. Similarly, we obtain the same upper bound for $E(B^{2})$. We finally conclude that

[TABLE]

As argued above, Theorem 13.5 in Billingsley (1999) now implies conditional tightness in probability given $\chi$.

The same holds true for $W_{\mathcal{I}}^{*}$ which, as was the case for the original estimator, is more easily shown than for $W_{{\mathcal{J}\mathcal{I}}}^{*}$ and it is thus left to the reader. All in all, conditional convergence in distribution of $\sqrt{n}(P_{\mathcal{I}\mathcal{J}}^{*}(s,\cdot)-\widehat{P}_{\mathcal{I}\mathcal{J}}(s,\cdot))$ on $D[s,\tau]$ given $\chi$ in probability follows by the functional delta-method. The Gaussian limit process is the same as for the original estimator.

Appendix E Proof of Theorem 3

Consider subsequences of increasing sample sizes for which $n_{1}/(n_{1}+n_{2})\rightarrow\lambda\in(0,1)$. Due to the independence of the eventual limit distribution on the specific value of $\lambda$, the convergence in distribution still holds under the bounded $\liminf$-$\limsup$-condition. The proof follows from a two-fold application of the continuous mapping theorem to the integration and the subtraction functional as well as from the independence of both samples. If $U$ is a zero-mean Gaussian process with covariance function $\Gamma_{{\mathcal{I}\mathcal{J}}}$, then $\int_{s}^{\tau}U(v)\mathrm{d}v$ has a zero-mean normal distribution with variance $\int_{s}^{\tau}\int_{s}^{\tau}\Gamma_{\mathcal{I}\mathcal{J}}(u,v)\mathrm{d}u\mathrm{d}v.$

Appendix F Detailed Derivation of the Covariance Function $\Gamma_{\mathcal{I}\mathcal{J}}$

As stated in Theorem 1, the covariance function $\Gamma_{\mathcal{I}\mathcal{J}}$ consists of several simpler covariance functions of which $cov(L_{0}(r),G_{\mathcal{J}|\mathcal{I}}(t))$ and $cov(L_{1}(r),G_{\mathcal{J}|\mathcal{I}}(t))$ still need to be determined. Note that, if $L_{2}$ denotes the limit Gaussian process of the normalized Aalen–Johansen estimator for the cumulative incidence function $F_{2}$ of the second risk, then

[TABLE]

due to $F_{0}+F_{1}+F_{2}=1$ on $[s,\tau]$. It is therefore enough to focus on the covariance function which involves the Aalen–Johansen limit. This may be derived based on the asymptotic linear representation of all three separate estimators as Hadamard-derivatives $\phi^{\prime}$, $\varphi^{\prime}$, and $\psi^{\prime}$, which are continuous and linear functionals, applied to an empirical process $\mathbb{P}_{n}$ and to $(\widehat{p}_{{\mathcal{J}\mathcal{I}}},\widehat{p}_{\mathcal{I}})$, respectively. More precisely, we have

[TABLE]

Here $\mathbb{P}_{n}$ is the empirical process based on the censoring or event times $T_{\ell}=\inf\{t\geq s:Y_{\ell}(t)=0\}$, the censoring or competing risk indicator $\delta_{\ell}=1\{T_{\ell 1}<T_{\ell 2}\}+2\cdot 1\{T_{\ell 1}>T_{\ell 2}\}$, and $\varepsilon_{\ell}=1\{X_{\ell}(s)\in\mathcal{I}\}$ as well as $\eta_{\ell}=1\{Y_{\ell}(s)=1\}$, where $T_{\ell 1}=\inf\{t\geq s:Z_{\ell}(t)=1,Y_{\ell}(t)=1\}$ and $T_{\ell 2}=\inf\{t\geq s:Z_{\ell}(t)=2,Y_{\ell}(t)=1\}$. Finally, $\mathbb{P}=P^{(T_{\ell},\delta_{\ell},\varepsilon_{\ell},\eta_{\ell})}$ is the distribution of each single observed quadruple. The empirical process is indexed by $f_{1z}(t,d,x,y)=1\{t\leq z,d=1,x=1,y=1\}$, $f_{2z}(t,d,x,y)=1\{t\leq z,d=2,x=1,y=1\}$, and $g_{z}(t,d,x,y)=1\{t\geq z,x=1,y=1\}$, $z\geq s$. This certainly is a Donsker class; see Example 3.9.19 in van der Vaart and Wellner (1996) for a similar indexing leading to the (all-cause) Nelson–Aalen estimator. In the following part, we calculate the exact structure of the above linear expansions for the case $n=1$. This suffices to determine the asymptotic covariance structure of the estimators as the covariances of the above linear functionals are the same for all $n$.

As argued in the main article, it does not matter whether the left- or right-continuous versions of the at-risk indicators $Y_{\ell}$ are considered. This is also true for the at risk functions involved in the Kaplan–Meier and Aalen–Johansen estimators. For technical convenience, we choose the right-continuous at-risk functions. Let $p_{\mathcal{A}_{\mathcal{J}},\mathcal{R}_{\mathcal{J}}}(t)=P(T_{1}\leq t,\delta_{1}\in\{1,2\},X_{1}(s)\in\mathcal{I},Y_{1}(s)=1)$, $p_{\mathcal{A}_{\mathcal{J}}}(t)=P(T_{1}\leq t,\delta_{1}=1,X_{1}(s)\in\mathcal{I},Y_{1}(s)=1)$, and $p_{\mathcal{R}_{\mathcal{J}}}(t)=P(T_{1}\leq t,\delta_{1}=2,X_{1}(s)\in\mathcal{I},Y_{1}(s)=1)$. The above Hadamard-derivatives are given by

[TABLE]

Thereof, $\phi_{1,-A_{0}}^{\prime}$ is the Hadamard-derivative of the product-integral leading to the Kaplan–Meier estimator, i.e.

[TABLE]

cf. e.g. Lemma 3.9.30 in van der Vaart and Wellner (1996), but take notice of the additional minus sign. Here, $A_{0}$ is the continuous cumulative (all-cause) hazard function for all individuals for which $X(s)\in\mathcal{I}$. The functional $\phi^{\prime}_{2,p_{\mathcal{A}_{\mathcal{J}},\mathcal{R}_{\mathcal{J}}},1/p_{\mathcal{I}}}$ is the Hadamard-derivative of the Wilcoxon functional, i.e.

[TABLE]

cf. e.g. Example 3.9.19 in van der Vaart and Wellner (1996). Finally, $\phi^{\prime}_{3,p_{\mathcal{A}_{\mathcal{J}},\mathcal{R}_{\mathcal{J}}},p_{\mathcal{I}}}$ is the Hadamard-derivative of $(a,b)\mapsto(a,1/b)$, after indexing the distribution with $(f_{1\cdot}+f_{2\cdot},g_{\cdot})$, i.e.

[TABLE]

From the above derivations, we have that

[TABLE]

We proceed similarly for the Aalen–Johansen estimator

[TABLE]

of $F_{1}(t)=\int_{s}^{t}F_{0}(u-)\mathrm{d}A_{1}(u)$, where now $A_{1}$ denotes the first cause-specific cumulative hazard of those individuals for which $X(s)\in\mathcal{I}$. Define $p_{\mathcal{A}_{\mathcal{J}}}(t)=P(T\leq t,\delta=1,X(s)\in\mathcal{I},Y(s)=1)$. Thus, as the estimator $\widehat{F}_{1}$ consists of the Wilcoxon functional applied to the left-continuous version of a product integral and another Wilcoxon functional (the Nelson–Aalen estimator), we have

[TABLE]

The tilde on top of the first Hadamard-derivative shall express the slight variation of the original functional, as we now use the left-continuous version of the integrand. Furthermore, $\tilde{\phi}^{\prime}_{3,p_{\mathcal{A}_{\mathcal{J}}},p_{\mathcal{I}}}(\mathbb{P}_{1}-\mathbb{P})(t)=(\mathbb{P}f_{1t}-p_{\mathcal{A}_{\mathcal{J}}}(t),-\frac{\mathbb{P}g_{t}-p_{\mathcal{I}}(t)}{p_{\mathcal{I}}^{2}(t)})$. Hence, the previous display can be expanded to

[TABLE]

Finally, we calculate the remaining linear functional applied to a pair of centered indicator functions that appear in the estimator $\widehat{p}_{\mathcal{J}|\mathcal{I}}$:

[TABLE]

For later calculations we note that

[TABLE]

that $P^{(T,\varepsilon,\eta)}(\mathrm{d}t,1,1)=-p_{\mathcal{I}}(\mathrm{d}t)$, and that

[TABLE]

where $P_{\mathcal{I}\mathcal{J}1}(s,t_{3},t)=P(X(t)=1,X(t_{3})\in\mathcal{J}\mid X(s)\in\mathcal{I})$.

After this preparation, we are now able to determine the asymptotic covariance. Because of the independencies and the identical distribution structure in the asymptotic, linear representation,

[TABLE]

We begin by considering the case $t_{2}\leq t_{3}$ for which the covariance reduces to

[TABLE]

The round bracket in the last line equals $p_{\mathcal{J}\mathcal{I}}(t_{3})-p_{\mathcal{J}|\mathcal{I}}(t_{3})p_{\mathcal{I}}(t_{3})=0$. Hence, asymptotic independence follows in case of $t_{2}\leq t_{3}$.

We hence continue with the case $t_{2}>t_{3}$. Let $f$ and $g$ be any suitable functions. The following side calculation due to a decomposition will prove beneficial:

[TABLE]

For a better presentation of the calculations, we split the covariance terms into two parts. First,

[TABLE]

Thereof, the sum of the last four terms can be summarized similarly as in the case $t_{2}\leq t_{3}$, i.e. it equals

[TABLE]

The sum of the remaining four terms equals, after integration by parts,

[TABLE]

Comparing these terms with those four which we derived first, we see that all integrals from $s$ to $t_{2}$ are cancelled out. Therefore, the following terms remain:

[TABLE]

Due to $P_{\mathcal{I}\mathcal{J}0}(s,t_{3},t_{3})=p_{{\mathcal{J}|\mathcal{I}}}(t_{3})F_{0}(t_{3})$, this equals

[TABLE]

We continue by deriving the second half of the asymptotic covariance function also in case of $t_{3}<t_{2}$, i.e.

[TABLE]

We apply integration by parts to the first integral. Note that, if $u$ is the variable of integration,

[TABLE]

Therefore, the first integral of $B$ equals

[TABLE]

All in all, the sum of both parts $A$ and $B$ equals

[TABLE]

which, due to $P_{\mathcal{I}\mathcal{J}1}(s,t_{3},t_{3})=F_{1}(t_{3})$, further simplifies to

[TABLE]

Appendix G Conditions (1)–(4) in the Simulation Study

In this section we prove that the conditions (1)–(4) of the main manuscript are satisfied in the situation of Section 4 therein. Recall that the stochastic development of the process $X$ depends on the state occupied at time $s=4$. Hence, conditional on $X(4)=1$ and after time $s=5$, $X$ is a Markov process with cumulative hazard matrix

[TABLE]

Otherwise, if $X(4)=0$, the process $X$ subsequently is a Markov process with cumulative hazard matrix

[TABLE]

Write $p=P(X(4)=1)$ and $q=1-p$. In this conditionally homogeneous Markov set-up it follows that the matrix $P(5,t)$ of transition probabilities at time $t\geq 5$ is given by

[TABLE]

where $\exp(D)=\sum_{k=0}^{\infty}\frac{1}{k!}D^{k}$ is the matrix exponential of a square matrix $D$. Now, it is obvious that $t\mapsto P(5,t)$ is infinitely often continuously differentiable, hence (1) and (2) are satisfied. In the same way, (3) and (4) are satisfied as well because $X$ is a Markov process conditionally on $X(4)$, so its state at an earlier intermediate time $t_{1}>s=5>4$ is irrelevant.