Nonparametric tests for transition probabilities in nonhomogeneous Markov processes

TL;DR

This paper introduces nonparametric two-sample tests for comparing transition probabilities in nonhomogeneous Markov processes, with rigorous significance assessment and extensions to complex cases, validated through simulations and real data application.

Contribution

It develops new nonparametric tests for transition probabilities in nonhomogeneous Markov processes, including extensions to complex scenarios and validation through simulations and real data.

Findings

Tests perform well with small samples

L2 and KS tests are consistent under fixed alternatives

Application to breast cancer treatment data demonstrates practical utility

Abstract

This paper proposes nonparametric two-sample tests for the direct comparison of the probabilities of a particular transition between states of a continuous time nonhomogeneous Markov process with a finite state space. The proposed tests are a linear nonparametric test, an L2-norm-based test and a Kolmogorov-Smirnov-type test. Significance level assessment is based on rigorous procedures, which are justified through the use of modern empirical process theory. Moreover, the L2-norm and the Kolmogorov-Smirnov-type tests are shown to be consistent for every fixed alternative hypothesis. The proposed tests are also extended to more complex situations such as cases with incompletely observed absorbing states and non-Markov processes. Simulation studies show that the test statistics perform well even with small sample sizes. Finally, the proposed tests are applied to data on the treatment of early breast cancer from the European Organization for Research and Treatment of Cancer (EORTC) trial 10854, under an illness-death model.

Figures

Tables

Table 1. Table 1: Simulation results about empirical type I error rates for the linear test (Linear), the L 2 superscript 𝐿 2 L^{2} -norm-based test ( L 2 superscript 𝐿 2 L^{2} ), and the Kolmogorov–Smirnov-type test (KS) under simulation scenaria 1 and 2.

			$\alpha =0.01$			$\alpha =0.05$
$H_0$ scenario	$n_1$	$n_2$	Linear	$L^2$	KS	Linear	$L^2$	KS
	50	50	0.009	0.007	0.015	0.051	0.054	0.063
	100	50	0.010	0.009	0.018	0.053	0.051	0.047
	100	100	0.012	0.011	0.012	0.060	0.060	0.054
	200	100	0.013	0.014	0.008	0.049	0.047	0.045
	200	200	0.010	0.012	0.010	0.051	0.047	0.052
	50	50	0.014	0.013	0.015	0.067	0.061	0.066
	100	50	0.014	0.014	0.017	0.057	0.053	0.065
	100	100	0.016	0.015	0.013	0.058	0.057	0.056
	200	100	0.011	0.014	0.009	0.058	0.060	0.057
	200	200	0.008	0.014	0.016	0.060	0.055	0.062

Table 2. Table 2: Simulation results about empirical type I error rates for the linear test (Linear), the L 2 superscript 𝐿 2 L^{2} -norm-based test ( L 2 superscript 𝐿 2 L^{2} ), and the Kolmogorov–Smirnov-type test (KS) under simulation scenaria 3 and 4.

			$\alpha =0.01$			$\alpha =0.05$
$H_0$ scenario	$n_1$	$n_2$	Linear	$L^2$	KS	Linear	$L^2$	KS
	50	50	0.016	0.014	0.011	0.075	0.069	0.066
	100	50	0.014	0.018	0.013	0.052	0.048	0.060
	100	100	0.011	0.012	0.012	0.048	0.047	0.046
	200	100	0.009	0.009	0.012	0.051	0.050	0.054
	200	200	0.013	0.015	0.013	0.053	0.051	0.058
	50	50	0.012	0.016	0.016	0.064	0.062	0.076
	100	50	0.012	0.012	0.013	0.060	0.059	0.067
	100	100	0.010	0.013	0.016	0.058	0.056	0.051
	200	100	0.007	0.008	0.010	0.044	0.054	0.058
	200	200	0.010	0.011	0.018	0.049	0.051	0.054

Table 3. Table 3: Simulation results about empirical power levels for the linear test (Linear), the L 2 superscript 𝐿 2 L^{2} -norm-based test ( L 2 superscript 𝐿 2 L^{2} ), and the Kolmogorov–Smirnov-type test (KS) under simulation scenaria 5 and 6.

			$\alpha =0.01$			$\alpha =0.05$
$H_1$ scenario	$n_1$	$n_2$	Linear	$L^2$	KS	Linear	$L^2$	KS
	50	50	0.074	0.069	0.055	0.170	0.162	0.166
	100	50	0.085	0.084	0.075	0.210	0.200	0.205
	100	100	0.106	0.101	0.101	0.251	0.240	0.240
	200	100	0.143	0.142	0.127	0.307	0.296	0.307
	200	200	0.233	0.206	0.207	0.448	0.417	0.397
	50	50	0.183	0.177	0.178	0.361	0.348	0.370
	100	50	0.246	0.230	0.228	0.442	0.438	0.445
	100	100	0.341	0.331	0.328	0.556	0.549	0.564
	200	100	0.458	0.464	0.480	0.703	0.696	0.708
	200	200	0.665	0.677	0.665	0.847	0.861	0.875

Table 4. Table 4: Simulation results about empirical rejection rates for the linear test (Linear), the L 2 superscript 𝐿 2 L^{2} distance test ( L 2 superscript 𝐿 2 L^{2} ), and the Kolmogorov–Smirnov-type test (KS) under simulation scenaria 7 and 8

			$\alpha =0.01$			$\alpha =0.05$
$H_1$ scenario	$n_1$	$n_2$	Linear	$L^2$	KS	Linear	$L^2$	KS
	50	50	0.015	0.033	0.033	0.064	0.099	0.123
	100	50	0.016	0.028	0.033	0.064	0.111	0.135
	100	100	0.017	0.056	0.060	0.074	0.164	0.169
	200	100	0.018	0.066	0.058	0.094	0.204	0.203
	200	200	0.026	0.102	0.105	0.095	0.304	0.288
	50	50	0.016	0.055	0.138	0.055	0.247	0.356
	100	50	0.018	0.091	0.189	0.064	0.326	0.441
	100	100	0.010	0.205	0.325	0.049	0.554	0.598
	200	100	0.016	0.351	0.490	0.062	0.712	0.795
	200	200	0.007	0.658	0.743	0.057	0.919	0.927

Equations

$$P_{hj}(s,t)$$

$${\hat{\mathbf{P}}}_n(s,t)={\text{Prodi}}_{(s,t]}\left[\mathbf{I}+d{\hat{\mathbf{A}}}_n(u)\right],s,t\in [0,\tau ],$$

$${\hat{A}}_{n,hj}(t)={\int }_0^t\frac{{\sum }_{i=1}^{n}d{N}_{ihj}(u)}{{\sum }_{i=1}^n{Y}_{ih}(u)},h\ne j.$$

$$D_{hj}(t)={\hat{W}}_{hj}(t)\left[{\hat{P}}_{n_1,hj}^{(1)}(0,t)-{\hat{P}}_{n_2,hj}^{(2)}(0,t)\right],t\in [0,\tau ]$$

$${\hat{W}}_{hj}(t)=\frac{{\bar{Y}}_h^{(1)}(t){\bar{Y}}_h^{(2)}(t)}{{\bar{Y}}_h^{(1)}(t)+{\bar{Y}}_h^{(2)}(t)}$$

$$Z_{hj}={\int }_{(0,\tau ]}{D}_{hj}(t)dm(t),$$

$$Z_{hj}={\int }_{(t_1,t_2]}{D}_{hj}(t)dm(t).$$

$$M_{ilm}^{(p)}(t)=N_{ilm}^{(p)}(t)-{\int }_{(0,t]}{Y}_{il}^{(p)}(u)d{A}_{0,lm}^{(p)}(u),$$

$$\sqrt{\frac{n_1{n}_2}{n_1+n_2}}{Z}_{hj}\stackrel{d}{\to }{G}_{hj},$$

$${\omega }_{hj}^2=(1-\lambda )E{\left[{\int }_{(0,\tau ]}{W}_{hj}(t){\gamma }_{1hj}^{(1)}(0,t)dm(t)\right]}^2+\lambda E{\left[{\int }_{(0,\tau ]}{W}_{hj}(t){\gamma }_{1hj}^{(2)}(0,t)dm(t)\right]}^2.$$

$${\gamma }_{ihj}^{(p)}(s,t)=\sum _{l\notin \mathcal{T}}\sum _{m\in \mathcal{I}}{\int }_s^t\frac{P_{0,hl}^{(p)}(s,u-)P_{0,mj}^{(p)}(u,t)}{E{Y}_{1l}^{(p)}(u)}d{M}_{ilm}^{(p)}(u),0\le s<t\le \tau ,p=1,2,$$

$${\hat{\omega }}_{hj}^2$$

$$\frac{Z_{hj}}{\hat{\omega}_{hj}\Big{/}\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}}.$$

$$Q_{1hj}={\left\{{\int }_{(0,\tau ]}{\left[D_{hj}(t)\right]}^{2}dm(t)\right\}}^{1/2}$$

$$Q_{2hj}=\underset{[0,\tau ]}{\sup }|D_{hj}(t)|.$$

$${\hat{B}}_{hj}(t)$$

$$\sqrt{\frac{n_1{n}_2}{n_1+n_2}}{D}_{hj}⇝\sqrt{1-\lambda }{\mathbb{G}}_{1hj}-\sqrt{\lambda }{\mathbb{G}}_{2hj},$$

$${\hat{B}}_{hj}⇝\sqrt{1-\lambda }{\mathbb{G}}_{1hj}-\sqrt{\lambda }{\mathbb{G}}_{2hj},$$

$${\sigma }_{hjp}(v,t)=E[W_{hj}(v){\gamma }_{1hj}^{(p)}(0,v)][W_{hj}(t){\gamma }_{1hj}^{(p)}(0,t)],p=1,2.$$

$$\sqrt{\frac{n_1{n}_2}{n_1+n_2}}{Q}_{1hj}\stackrel{d}{\to }{\left\{{\int }_{(0,\tau ]}{\left[\sqrt{1-\lambda }{\mathbb{G}}_{1hj}(t)-\sqrt{\lambda }{\mathbb{G}}_{2hj}(t)\right]}^{2}dm(t)\right\}}^{1/2},$$

$$\sqrt{\frac{n_1{n}_2}{n_1+n_2}}{Q}_{2hj}\stackrel{d}{\to }\underset{t\in [0,\tau ]}{\sup }\left|\sqrt{1-\lambda }{\mathbb{G}}_{1hj}(t)-\sqrt{\lambda }{\mathbb{G}}_{2hj}(t)\right|.$$

$${\left\{{\int }_{(0,\tau ]}{\left[{\hat{B}}_{hj,r}(t)\right]}^{2}dm(t)\right\}}^{1/2},r=1,\dots ,R$$

$${\hat{B}}_{hj,r}(t)$$

$${\tilde{Z}}_{hj}={\int }_H{\tilde{D}}_{hj}(t)dm(t),$$

$${\tilde{D}}_{hj}(t)={\hat{W}}_{hj}(t)\left[{\tilde{P}}_{n_1,hj}^{(1)}(s,t)-{\tilde{P}}_{n_2,hj}^{(2)}(s,t)\right],t\in H,$$

$${\tilde{Q}}_{1hj}={\left\{{\int }_H{\left[{\tilde{D}}_{hj}(t)\right]}^{2}dm(t)\right\}}^{1/2}$$

$${\tilde{Q}}_{2hj}=\underset{t\in H}{\sup }|{\tilde{D}}_{hj}(t)|.$$

$$\underset{t\in H}{\sup }|{\tilde{P}}_{n,hj}(s,t)-P_{hj}(s,t)|\stackrel{p}{\to }0,$$

$$\sqrt{n}[{\tilde{P}}_{n,hj}(s,t)-P_{hj}(s,t)]=\frac{1}{\sqrt{n}}\sum _{i=1}^n{\varphi }_{ihj}(s,t)+o_p(1),$$

$$\underset{t\in H}{\sup }\left|n^{-1/2}\sum _{i=1}^n[{\hat{\varphi }}_{ihj}(s,t)-{\varphi }_{ihj}(s,t)]{\xi }_i\right|\stackrel{p}{\to }0,$$

Full text

Keywords: Competing risks; Crossing curves; Illness-death model; Missing absorbing states; Multistate model.

1 Introduction

Continuous time nonhomogeneous Markov processes with a finite state space are important in many areas of science and particularly in medicine and public health (Tattar and Vaman, 2014; Bakoyannis et al., 2019). Consideration of specific transitions between two states of a multi-state process can provide a deeper and more detailed insight about the treatment effect in clincal trials compared to the analysis of standard survival outcomes, such as event-free survival (Le-Rademacher et al., 2018). Important special cases of a Markov process are the univariate survival model, the competing risks model, and the Markov illness-death model (Andersen et al., 2012).

The stochastic behaviour of a Markov process can be described by either the transition intensities, which represent the instantaneous rates of transition between two states, or the transition probabilities. The transition probabilities are also known as survival functions in the framework of the univariate survival model, and as cumulative incidence functions in the competing risks model. It is important to note that, in general, a difference in the transition intensities between two groups does not necessarily imply a difference in the corresponding transition probabilities and vice versa. This phenomenon has been well documented for the special case of the competing risks model (Gray, 1988; Pepe, 1991; Putter et al., 2007; Bakoyannis and Touloumi, 2012). Nonparametric tests for comparing transition intesities between groups in general Markov multi-state processes have been well developed (Andersen et al., 2012). However, the issue of nonparametric comparison of transition probabilities in general Markov multi-state processes has not received much attention. Nevertheless, transition probabilities, unlike transition intensities, directly quantify clinical prognosis (Bakoyannis et al., 2019), which is the target of scientific interest in many applications.

Nonparametric estimation of the transition probabilities of a general Markov process can be performed using the Aalen–Johansen estimator (Aalen and Johansen, 1978). The issue of nonparametric comparison of transition probabilities under the univariate survival model has be extensively studied in the literature. For a review of these methods see Kalbfleisch and Prentice (2011) and Andersen et al. (2012). A number of researchers have proposed nonparametric tests for the comparison of transition probabilities for the special case of the competing risks model (Gray, 1988; Pepe and Mori, 1993; Lin, 1997). Dabrowska and Ho (2000) proposed a graphical procedure based on simultaneous confidence bands to test for differences between transition probabilities in a general Markov process. However, their method imposes proportional hazards assumptions for the transition intensities and, thus, it is not fully nonparametric. Also, this approach does not provide the actual level of statistical significance. Tattar and Vaman (2014) proposed two nonparametric tests for the comparison of the whole transition probability matrices between $k$ groups, by comparing all the possible transition intensities. The first test only compares the transition probability matrices at a specific time point $t_{0}$, while the second test is a Kolmogorov–Smirnov-type test based on the supremum norm. However, the tests proposed by Tattar and Vaman (2014) do not provide a direct comparison of the transition probability of a particular transition, which is frequently of scientific interest (Le-Rademacher et al., 2018). A statistically significant result with these tests only indicates a difference in any transition between groups. Recently, Bluhmki et al. (2018) proposed a wild bootstrap approach for the Aalen–Johansen estimator, which can be used to construct a simultaneous confidence band for the difference between the transition probabilities of two independent groups. This approach, which is related to a Kolmogorov–Smirnov-type test, can be used as a graphical two-sample comparison procedure at a predetermined $\alpha$ level. However, this approach does not provide the actual level of statistical significance and, also, a Komogorov–Smirnov-type test may not be the most powerful nonparametric test for every alterantive hypotheses. Additionally, there is no rigorous justification about the consistency of this graphical hypothesis testing procedure against any fixed alternative hypothesis (Van der Vaart, 2000). Last but not least, the proposed approach is not readily adaptable to more complex situations such as cases with missing data.

This paper addresses the issue of direct nonparametric two-sample comparison of the transition probabilities of a particular transition in a general continuous-time nonhomogeneous Markov process with a finite state space. For this, we propose a linear nonparametric test, an $L^{2}$-norm-based test and a Kolmogorov–Smirnov-type test. The asymptotic null distributions of the tests are derived. The evaluation of the actual level of statistical significance is based on rigorous procedures justified through the use of modern empirical process theory. Moreover, the $L^{2}$-norm-based and Kolmogorov–Smirnov-type tests are shown to be consistent against any fixed alternative hypothesis (Van der Vaart, 2000). We also propose extensions related to interesting partical problems such as cases with missing absorbing states (Bakoyannis et al., 2019) and non-Markov processes (Putter and Spitoni, 2018). The proposed tests exhibit good small sample properties as illustrated in our simulation experiments. Finally, the tests are applied to data on the treatment of early breast cancer from the European Organization for Research and Treatment of Cancer (EORTC) trial 10854.

Compared to the previous work by Bluhmki et al. (2018), which used counting process theory arguments in their derivations, we justify the properties of the proposed tests through the use of modern empirical process theory (Van Der Vaart and Wellner, 1996; Kosorok, 2008). As it will be argued later in the text, the practical advantage of our derivations lies on the fact that our proposed tests can be straightforwardly adapted to more complex settings such as cases with incompletely observed absorbing states (Bakoyannis et al., 2019). This can be done by replacing the influence function of the standard Aalen–Johansen estimator with the influence function of any other well-behaved and asymptotically linear estimator of the transition probabilities in our proposed testing procedures. Such adaptations are not trivial within the framework of the graphical testing procedure proposed by Bluhmki et al. (2018). An important reason for this is that with more complex estimators, certain predictability conditions assumed by counting process and martingale theory techniques are violated. For such situations, empirical process theory provides a powerful alternative tool. Moreover, we provide two additional tests, a linear test and an $L^{2}$-norm-based test, which may be more powerful compared to a Kolmogorov–Smirnov-type test in certain settings. Additionally, we argue about the consistency of our $L^{2}$-norm-based and Kolmogorov–Smirnov-type tests against any fixed alternative hypothesis. Finally, our tests provide the actual level of statistical significance which is useful in pactical applications.

The structure of this paper is as follows. In Section 2 we introduce some notation about Markov processes, provide the proposed nonparametric tests, and consider extensions to more complex situations that are frequently met in practice. Section 3 presents a simulation study to evaluate the small sample performance of the proposed tests. Section 4 illustrates the use of the proposed tets using data from the EORTC trial 10854. Finally, Section 5 conlcudes the article with some key remarks. Outlines of the asymptotic theory proofs are provided in the Appendix.

2 Two-sample nonparametric tests

2.1 Nonparametric estimation of transition probabilities

The stochastic behaviour of a Markov process $\{X(t):t\geq 0\}$ with a finite state space $\mathcal{I}=\{1,\ldots,q\}$ can be described by the $q\times q$ transition probability matrix $\mathbf{P}_{0}(s,t)=(P_{hj}(s,t))$ whose elements are the transition probabilities

[TABLE]

where \mathcal{F}_{s^{-}}=\sigma\big{\langle}\{N_{hj}(u):0\leq u<s,h\neq j\}\big{\rangle} is the event history prior to time $s$, with $N_{hj}(t)$ being the number of direct transitions from state $h$ to state $j$, $h\neq j$, in $[0,t]$, $\tau=\sup\{t:\int_{0}^{t}a_{hj}(u)du\equiv A_{hj}(t)<\infty,h\neq j\}$, and $a_{hj}(t)=\lim_{h\downarrow 0}P_{hj}(t,t+h)/h$, $h\neq j$, is the transition intensity at time $t$. The conditional independence between the probability of $X(t)$ and the prior history $\mathcal{F}_{s^{-}}$, conditionally on $X(s)$, is the so-called Markov assumption. Because $\mathbf{P}_{0}(s,t)$ is a stochastic matrix we have that $a_{hh}(t)=-\sum_{j\neq h}a_{hj}(t)$.

The observed data from a sample of i.i.d. observations of a Markov process are the counting processes $\{N_{ihj}(t):h\neq j,t\in[0,\tau]\}$, $i=1,\ldots,n$, which represent the number of direct transitions of the $i$th observation from the state $h$ to the state $j$ by time $t$, and the at-risk processes $\{Y_{ih}(t):h\in\mathcal{I},t\in[0,\tau]\}$ which are the indicator functions of whether the $i$th observation is at the state $h\in\mathcal{I}$ just before time $t\in[0,\tau]$. Based on such a sample, the transition probability matrix of a nonhomogeneous Markov process can be estimated using the Aalen–Johansen estimator (Aalen and Johansen, 1978):

[TABLE]

where $\prodi$ is the product integral and $\hat{\mathbf{A}}_{n}(t)$ a $q\times q$ matrix whose elements are the Nelson–Aalen estimates of the cumulative transition intensities

[TABLE]

2.2 Linear nonparametric tests

First consider the two-sample problem of comparing the transition probabilities $P_{0,hj}^{(1)}(s,\cdot)$ and $P_{0,hj}^{(2)}(s,\cdot)$, $s\in[0,\tau)$, of two populations of interest, for a particular transition $h\rightarrow j$, with $h,j\in\mathcal{I}$. For simplicity of presentation we will set the starting point $s=0$ for the remainder of the paper. Based on two independent random samples of $n_{1}$ and $n_{2}$ observations from the two populations, define the pointwise weighted difference

[TABLE]

where $\hat{W}_{hj}(t)$ is a weight function and $\hat{P}_{n_{1},hj}^{(1)}(0,t)$ and $\hat{P}_{n_{2},hj}^{(2)}(0,t)$ are the nonparametric Aalen–Johansen estimates of the transition probabilities of the two populations under comparison. Example of weight function choices are $W_{hj}(t)=1$ and

[TABLE]

where $\bar{Y}_{h}^{(p)}(t)=n_{p}^{-1}\sum_{i=1}^{n_{p}}Y_{ih}^{(p)}(t)$, $p=1,2$. The latter choice assigns more weight to times with more observations at risk. A natural linear test for the null hypothesis $H_{0}:P_{0,hj}^{(1)}=P_{0,hj}^{(2)}$ is the area under the weighted difference curve

[TABLE]

where $m$ is the Lebesgue measure on the Borel $\sigma$-algebra $\mathcal{B}([0,\tau])$. To establish the asymptotic distribution of the test statistic $Z_{hj}$ we assume the following conditions.

• C1.

The potential right censoring and left truncation are independent of the counting processes $\{N_{hj}(t):h\neq j,t\in[0,\tau]\}$ and noninformative about $\mathbf{P}_{0}(s,t)$.

• C2.

$n_{1}/(n_{1}+n_{2})\rightarrow\lambda\in(0,1)$ as $\min(n_{1},n_{2})\rightarrow\infty$.

• C3.

The counting processes $\{N_{hj}(t):h\neq j,t\in[0,\tau]\}$ are bounded in the sense that $\Pr(N_{hj}(\tau)\leq C)=1$ for some constant $C\in(0,\infty)$.

• C4.

$\inf_{t\in[0,\tau]}E[Y_{h}(t)]>0$ for all the transient states $h\in\mathcal{I}$.

• C5.

The cumulative transition intensities $\{A_{0,hj}(t):h\neq j,t\in[0,\tau]\}$ are continuous functions of bounded variation on $[0,\tau]$.

• C6.

The weight $\hat{W}_{hj}(t)$ converges uniformly to a nonnegative uniformly bounded function $W_{hj}(t)$ on $[0,\tau]$.

Remark 1**.**

In some applications condition C4 may not be satisfied for some timepoints for one or more states $h\in\mathcal{I}$. In such cases one can restrict the comparison interval to $[t_{1},t_{2}]$ with $0<t_{1}<t_{2}<\tau$, such that $\inf_{t\in[t_{1},t_{2}]}E[Y_{h}(t)]>0$ for those $h\in\mathcal{I}$. In such cases the test statistic becomes

[TABLE]

Before stating the theorem about the asymptotic distribution of test statistic we define the functions

[TABLE]

where $N_{ilm}^{(p)}(t)$ and $Y_{il}^{(p)}(t)$ are the counting and at-risk processes of the $i$th observation in the $p$th sample at time $t$. Also, define $\mathcal{T}$ to be the subset of $\mathcal{I}$ which contains the potential absorbing states. The set $\mathcal{T}$ will be null for non-absorbing Markov processes.

Theorem 1 provides the asymptotic distribution of $Z_{hj}$ under the null hypothesis.

Theorem 1**.**

Suppose that conditions C1-C6 hold. Then under the null hypothesis

[TABLE]

where $G_{hj}\sim N(0,\omega_{hj}^{2})$ and

[TABLE]

with

[TABLE]

for $i=1,\ldots,n_{p}$.

Remark 2**.**

The functions $\gamma_{ihj}^{(p)}(s,t)$, p=1,2, in Theorem 1 are the influence functions of the Aalen–Johansen estimator.

A consistent (in probability) estimator of the variance $\omega_{hj}^{2}$ is

[TABLE]

where $\gamma_{ihj}^{(p)}(0,t)$, $p=1,2$, are estimated by replacing the expectations with sample averages and the unknown parameters with their uniform consistent estimates. Now, Theorem 1 and $\hat{\omega}_{hj}$ can be used to constuct a $Z$-test for the null hypothesis as:

[TABLE]

The actual significance level can then be avaluated under the standard normal distribution as usual.

2.3 $L^{2}$-norm-based and Kolmogorov–Smirnov-type tests

A linear test is not the optimal choice when the two transition probability curves under comparison cross at one or more timepoints. In this section, we propose alternative tests for such situations. The first test is a test based on an $L^{2}$ norm

[TABLE]

while the second test is a Kolmogorov–Smirnov-type test

[TABLE]

The Kolmogorov–Smirnov-type test is related to the graphical hypothesis testing procedure proposed by Bluhmki et al. (2018). The asymptotic null distributions of these tests are complicated. However, significance level can be easily calculated numerically by proper simulation realizations from the null distribution of these test statistics. Theorem 2 provides the basis for an approach to properly simulate realizations from the null distributions of $Q_{1hj}$ and $Q_{2hj}$. Before stating Theorem 2 define the estimated functions

[TABLE]

where $\{\xi_{i}^{(p)}\}_{i=1}^{n_{p}}$, $p=1,2$, are independent draws from $N(0,1)$.

Theorem 2**.**

Suppose that conditions C1-C6 hold. Then under the null hypothesis

[TABLE]

and, conditionally on the observed data,

[TABLE]

where $\mathbb{G}_{1hj}$ and $\mathbb{G}_{2hj}$ are two independent tight zero-mean Gaussian processes with covariance functions

[TABLE]

Corollary 1**.**

By Theorem 2 and the continuous mapping theorem it follows that under the null hypothesis

[TABLE]

and

[TABLE]

The asymptotic null distributions of the omnibus tests are quite complicated and, thus, they are of limited use in terms of evaluating the significance level. However, Theorem 2 provides justification about a way to numerically calculate $p$-values through a simple simulation technique. This can be performed as follows. In light of Theorem 2, one can simulate from the asymptotic null asymptotic distributions of the tests $Q_{1hj}$ and $Q_{2hj}$ by simulating multiple versions of $\{\xi_{ir}^{(1)}\}_{i=1}^{n_{1}}$ and $\{\xi_{ir}^{(2)}\}_{i=1}^{n_{2}}$ independently from $N(0,1)$ for $r=1,\ldots,R$, and then calculating a sample for the above null distributions as

[TABLE]

and $\sup_{t\in[0,\tau]}\left|\hat{B}_{hj,r}(t)\right|$, $r=1,\ldots,R$, respectively, where

[TABLE]

Now, the significance level for each test can be calculated as the proportion of realizations from the corresponding null distribution that is greater than or equal to the calculated tests statistic value from the observed data.

The tests $Q_{1hj}$ and $Q_{2hj}$ are consistent for every fixed alternative hypothesis with $P_{0,hj}^{(1)}\neq P_{0,hj}^{(2)}$. This follows from Theorem 2, the uniform consistency of the Aalen–Johansen estimator of the transition probabilities (Aalen and Johansen, 1978), condition C6, the continuity of these tests in $D_{hj}(t)$, and Lemma 14.15 in Van der Vaart (2000).

2.4 Extensions to more complex settings

Many complications that frequently occur in practice make the application of the proposed tests improper. An important example is the problem of incompletely observed absorbing states, where missingness occurs either due to the usual nonresponse or the study design (Bakoyannis et al., 2019). A special case of this is the issue of missing causes of death in biomedical applications. In such cases, a complete case analysis, which discards cases with a missing cause of death, is well known to lead to biased estimates (Gao and Tsiatis, 2005; Lu and Liang, 2008; Bakoyannis et al., 2019). In general, more complicated cases require extensions of the standard Aalen–Johansen estimator, denoted by $\tilde{P}_{n,hj}(s,t)$, to consistently estimate the transition probabilities of interest over a compact interval $H\subset[0,\tau]$. In such cases, one can replace the standard Aalen–Johansen estimator with another approriate estimator $\tilde{P}_{n,hj}(s,t)$ in the testing procedures. Then, the linear test becomes

[TABLE]

where

[TABLE]

while the $L^{2}$-norm based and Kolmogorov–Smirnov-type tests become

[TABLE]

and

[TABLE]

The following conditions ensure the validity of the proposed testing procedures in more complex settings.

• D1.

The estimator $\tilde{P}_{n,hj}(s,\cdot)$ is consistent in the sense

[TABLE]

for some $s\geq 0$, where $H$ is a compact subset of $[0,\tau]$.

• D2.

The estimator $\tilde{P}_{n,hj}(s,\cdot)$ is an asymptotically linear estimator with

[TABLE]

where the influence functions $\phi_{ihj}(s,t)$ belong to a Donsker class indexed by $H$.

• D3.

The empirical versions of the influence functions $\hat{\phi}_{ihj}(s,t)$ satisfy

[TABLE]

where $\xi_{i}$ are independent random draws from $N(0,1)$.

Remark 3**.**

Condition D2 is sufficient for establishing the weak convergence of the estimator $\tilde{P}_{n,hj}(s,\cdot)$ to a tight mean-zero Gaussian process. Condition D3 along with the conditional multiplier central limit theorem (Van Der Vaart and Wellner, 1996; Kosorok, 2008) and condition D2, provide a simulation approach for the construction of simultaneous confidence bands (Kosorok, 2008). Therefore, conditions D1-D3 are expected to have been established in works extending the standard Aalen–Johansen estimator to more complex settings. This is the case, for example, for the nonparametric estimator of the transition probability matrix with incompletely observed absorbing states (Bakoyannis et al., 2019).

Hypothesis testing in more complex settings can be simply performed by replacing the influence functions $\gamma_{ihj}^{(p)}(s,t)$, $p=1,2$, of the standard Aalen–Johansen estimator with the influence functions $\phi_{ihj}^{(p)}(s,t)$ of the estimator $\tilde{P}_{n,hj}(s,t)$. The theorems stated below justify the direct use of the proposed tests in more complex situations. Before stating those theorems define the functions

[TABLE]

where $\{\xi_{i}^{(p)}\}_{i=1}^{n_{p}}$, $p=1,2$, are independent draws from $N(0,1)$.

Theorem 3**.**

Suppose that conditions C2, C6, D1 and D2 hold. Then under the null hypothesis

[TABLE]

where $\tilde{G}_{hj}\sim N(0,\theta_{hj}^{2})$ and

[TABLE]

The proof of Theorem 3 involves the same arguments to those used in the proof of Theorem 1 given in the Appendix.

Theorem 4**.**

Assume that conditions C2, C6, and D1–D3 are satisfied. Then, under the null hypothesis

[TABLE]

and, conditionally on the observed data,

[TABLE]

where $\tilde{\mathbb{G}}_{1hj}$ and $\tilde{\mathbb{G}}_{2hj}$ are two independent tight zero-mean Gaussian processes with covariance functions

[TABLE]

The proof of Theorem 4 follows from similar arguments to those used in the proof of Theorem 2 given in the Appendix.

2.4.1 Missing absorbing states

In many settings one can observe that a process has arrived at some absorbing state, but the actual absorbing state is unobserved for some study participants, such as in cases with missing causes of death. For such situations, Bakoyannis et al. (2019) proposed a nonparametric maximum pseudolikelihood estimator (NPMPLE) under a missing at random assumption. To review this estimator, let $\Delta_{ij}$ be an indicator variable with $\Delta_{ij}=1$ if the $i$th observation arrived at the absorbing state $j\in\mathcal{T}$, and $\Delta_{ij}=0$ otherwise. Also, let $R_{i}$ be another indicator variable with $R_{i}=1$ indicating that the absorbing state of the $i$th observation has been successfully ascertained. Finally, let $\pi_{j}(\mathbf{O}_{i},\boldsymbol{\beta}_{0})$ be the probability that $\Delta_{ij}=1$ given the fully observed data $\mathbf{O}_{i}$, under a parametric model indexed by an unknown Euclidean parameter $\boldsymbol{\beta}_{0}$. In this setting, the cumulative transition intensities can be estimated using the NPMPLE:

[TABLE]

where

[TABLE]

with $\hat{\boldsymbol{\beta}}_{n}$ being a consistent estimator of $\boldsymbol{\beta}_{0}$. The transition probability matrix can then be estimated as

[TABLE]

where the components of the matrix $\tilde{\mathbf{A}}_{n}(u)$ are $\tilde{A}_{n,hj}(u)$. By Theorems 1 and 2 in Bakoyannis et al. (2019) and calculations provided in the proof of Theorem 2 in the same source, the NPMPLE estimator satisfies the conditions D1-D3 above. Therefore, if the conditions in Bakoyannis et al. (2019) and the conditions C2 and C6 above are satisfied, two-sample comparison can be performed by utilizing the NPMPLE of the transition probabilities along with the corresponding influence functions in the proposed tests. This is justified by Theorems 3 and 4 above.

2.4.2 Non-Markov processes

Trivially, the Aalen–Johansen estimator $\hat{P}_{n,hj}(0,\cdot)$ is uniformly consistent for the transition probability $P_{hj}(0,\cdot)$ even under a non-Markov process (Datta and Satten, 2001; Titman, 2015). When the interest lies on the marginal $\Pr(X(t)=j|X(s)=h)$, i.e. unconditionally on the prior history $\mathcal{F}_{s^{-}}$, for some $s>0$, under a non-Markov process, then the landmark Aalen–Johansen estimator is consistent for $\Pr(X(t)=j|X(s)=h)$ (Putter and Spitoni, 2018) under the conditions of Datta and Satten (2001) and, also, the assumption that $\Pr(X(s)=h)>0$. The landmark Aalen–Johansen estimator is essentially equivalent to the standard Aalen–Johansen estimator, except for the fact that only observations with $X(s)=h$ are considered. This is achieved by considering the modified counting and at-risk processes $\tilde{N}_{ihj}(t)=N_{ihj}(t)I(X(s)=h)$ and $\tilde{Y}_{ih}(t)=Y_{ih}(t)I(X(s)=h)$, for $t\geq s$. Therefore, the influence functions of the landmark Aalen–Johansen estimator are the same to that of the standard Aalen–Johansen estimator, with the only exception that the former involves the modified $\tilde{N}_{ihj}(t)$ and $\tilde{Y}_{ih}(t)$ instead of the standard counting and at-risk processes $N_{ihj}(t)$ and $Y_{ih}(t)$. Consequently, it is clear that conditions D1–D3 are satified if $\Pr(X(s)=h)>0$ and, also, if the conditions in Datta and Satten (2001) hold. Thus, in light of Theorems 3 and 4, the proposed nonparametric tests can be used with non-Markov processes by utilizing the landmark Aalen–Johansen estimator.

3 Simulation studies

To evaluate the finite sample performance of the proposed test statistics, we conducted a simulation study. We considered a nonhomogeneous Markov process with 2 transient states $\{1,2\}$ and 1 absorbing state $\{3\}$, under the illness-death model without recovery (Andersen et al., 2012). This model is illustrated in Figure 1. In this simulation study, we focused on the null hypothesis $H_{0}:P_{12}^{(1)}=P_{12}^{(2)}$. Initially, we independently generated the times from state 1 to states 2 and 3 by assuming the cumulative transition intensities $A_{12}(t)=\alpha_{1}t$ and $A_{13}(t)=t/2$. For observations that first arrived at the transient state 2, we generated the time from state 2 to the absorbing state 3, assuming a cumulative transition intensity $A_{23}(t)=\alpha_{2}t$. The parameter values considered were $\alpha_{1}\in\{0.4,0.6,1.2,1.4\}$ and $\alpha_{2}\in\{0.25,0.5,0.75\}$. Then, the right censoring times were independently simulated from Exp(0.25). Under this set-up the transition probability of interest was

[TABLE]

Different simulation scenarios were considered according to the sample sizes $n_{p}$, $p=1,2,$ and the parameter values of the two groups. 1,000 datasets were simulated for each scenario, and the $L^{2}$ distance test and Kolmogorov–Smirnov-type test were calculated using 1,000 independent simulations of $\{\xi_{i}^{(1)}\}_{i=1}^{n_{1}}$ and $\{\xi_{i}^{(2)}\}_{i=1}^{n_{2}}$ from $N(0,1)$. Finally, the weight function $W_{hj}(t)=1$ was considered in all cases.

Simulation results regarding the empirical type I error rates are presented in Tables 1 and 2, respectively. Under these scenaria, the empirical type I errors rates for all tests were close to the nominal $\alpha$ levels, even in situations with small sample sizes. Thus, these results provide numerical evidence for the validity of the proposed hypothesis testing procedures under $H_{0}$. Simulation results regarding the empirical power levels under alternative hypotheses with non-crossing transition probabilities are presented in Table 3. Under these scenaria, the empirical power levels increased with sample size and, also, with a more pronounced difference between the two groups, as expected. The power levels for the three tests were in general similar. However, under a less pronounced difference between the two groups, the linear test exhibited a somewhat larger empirical power with larger sample sizes. These results provide numerical evidence for the consistency of the proposed tests with non-crossing transition probabilities. Simulation results regarding the empirical power levels under alternative hypotheses with crossing transition probabilities are presented in Table 4. These scenaria illustrate numerically the inconsistency of the linear tests with crossing transition probabilities, as the empirical power levels did not systematically increase with sample size. On the contrary, the empirical power of the $L^{2}$-norm-based and Kolmogorov–Smirnov-type tests increased with sample size and with a more pronounced difference between the two groups. These results indicate numerically the consistency of the omnibus tests against alternatives with crossing transition probabilities.

4 Data analysis

In this section we analyze the data on treatment of early breast cancer from the European Organization for Research and Treatment of Cancer (EORTC) trial 10854. This randomized clinical trial was conducted to evaluate whether the combination of surgery with polychemotherapy is benefical to early breast cancer patients compared to surgery alone. The original analysis of this clinical trial was presented in Van der Hage et al. (2001).

In this trial, 1619 patients where randomly assinged to the surgery group and 1559 to the surgery plus polychemotherapy group. The data set contains information about the time to cancer relapse or death. Therefore, an illness-death model is a natural choice for this data set. It is important to note that the transition probability to relapse, which was not analyzed in the original analysis of this trial, is a non-monotonic function of time as patients can move to the “death” state after relapse. Thus, standard survival and competing risks analysis methods are not applicable for this transition probability. Here, we focus on this probability which can be interpreted as the probability of being alive and in relapse. The estimated transition probabilities of relapse in the two intervention groups are presented in Figure 2. Based on Figure 2, the probability of being alive and in relapse was lower in the group that received polychemotherapy during surgery. To perform hypothesis testing here we considered the weight function $W_{12}(t)=1$. For the $L^{2}$-norm-based and Kolmogorov–Smirnov-type tests we considered 1,000 standard normal simulation realizations. The $p$-value from the linear test was 0.001, while the $p$-values from the $L^{2}$-norm-based and Kolmogorov–Smirnov-type tests were $<$0.001 and 0.004, respectively. These results provide evidence for the superiority of the surgery plus polychemotherapy combination with respect to the transition probability of relapse, in early breast cancer patients.

5 Concluding remarks

This paper addressed the issue of direct nonparametric two-sample comparison of transition probabilities $P_{hj}(s,\cdot)$, $s\in[0,\tau)$, for a particular transition $h\rightarrow j$ in a continuous time nonhomogeneous Markov process with a finite state space. The proposed tests were a linear nonparametric test, an $L^{2}$-norm-based test and a Kolmogorov–Smirnov-type test. Rigorous approaches to evaluate the significance level grounded on modern empirical process theory were provided. Moreover, the $L^{2}$-norm-based and Kolmogorov–Smirnov-type tests were argued to be consistent against any fixed alternative hypothesis. We also considered extensions of the tests to more complex situations such as cases with missing absorbing states (Bakoyannis et al., 2019) and non-Markov processes (Putter and Spitoni, 2018). The simulation study provided numerical evidence for the validity of the proposed testing procedures, which exhibited good performance even with small sample sizes. Finally, a data analysis of a clinical trial on early breast cancer illustrated the utility of the proposed tests in practice.

The issue of nonparametric comparison of transition probabilities in general nonhomogeneous Markov processes has received little attention in the literature. To the best of our knowledge, the only fully nonparametric approach for comparing the transitions probabilitis for a particular transtion in general non-homogeneous Markov processes is a graphical procedure proposed by Bluhmki et al. (2018). This proposal is based on the construction of a simultaneous confidence band for the difference between the transition probabilities of two groups. However, this approach does not provide the exact level of statistical significance. Also, the justification of this approach was based on counting process theory arguments and not on modern empirical process theory. A concequence of that is that this approach cannot be directly adapted to more complex settings that are frequently occur in practice, such as cases with missing absorbing states. An important reason for this is that with more complex estimators, certain predictability conditions assumed by counting process and martingale theory techniques are violated. On the contrary, our proposed methods can be trivially adapted to many other complex settings, provided that appropriate estimators, in the sense of conditions D1–D3, of the transition probabilities exist. Such adaptations can be theoretically justified using the Theorems 3 and 4 provided in our manuscript.

The proposed tests can be easily adapted for the comparison of state occupation probabilities $\Pr(X(t)=j)\equiv P_{j}(t)=\sum_{h\in\mathcal{I}}P_{h}(0)P_{hj}(0,t)$, as these are simple linear combinations of the transition probabilities. The state occupation probabilities describe the marginal behavior, i.e. unconditional on the prior history, of the processes and are of interest in many applications, such as in HIV studies focusing on the event history of patients in HIV care (Lee et al., 2018). It is important to note that these probabilities can be consistently estimated using the Aalen–Johansen estimator even in non-Markov processes (Datta and Satten, 2001). It is not hard to justify conditions D1–D3 for the state occupation probabiltities under a set of weak regularity conditions. Thus, Theorems 3 and 4 provide a rigorous justification about the use of the proposed tests for comparing state occupation probabilities.

Acknowledgement

This project was supported, in part, by the Indiana Clinical and Translational Sciences Institute funded, in part by Grant Number UL1TR002529 from the National Institutes of Health, National Center for Advancing Translational Sciences, Clinical and Translational Sciences Award. We would like to thank the European Organisation for Research and Treatment of Cancer (EORTC) for sharing with us the data from the EORTC trial 10854. The content of this manuscript is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health and the EORTC.

Appendix A Outlines of proofs

Outlines of the proofs of Theorems 1 and 2 are provided below. The proofs of Theorems 3 and 4 follow from similar arguments and, therefore, are omitted. The proofs rely on empirical process theory techniques (Van Der Vaart and Wellner, 1996; Kosorok, 2008). Before providing the proofs it is useful to introduce some notation. First, let $\mathcal{O}$ be the sample space, and $O$ an arbitraty sample point in $\mathcal{O}$. Now, define $\mathbb{P}_{n}f=\frac{1}{n}\sum_{i=1}^{n}f(O_{i})$, for some measurable function $f:\mathcal{O}\mapsto\mathbb{R}$. Also, define $Pf=\int_{\mathcal{O}}fdP$ to be the expectation of $f$ under the probability measure $P$ on the measurable space $(\mathcal{O},\mathcal{A})$, where $\mathcal{A}$ is a $\sigma$-algebra on $\mathcal{O}$. For simplicity, but without loss of generality, we set the starting point $s=0$ in the following proofs. It has to be noted that conditions C1 and C3–C5 imply the uniform consistency of the standard Aalen–Johansen estimator. This can be shown using similar arguments to those used in the proof of Theorem 1 in Bakoyannis et al. (2019).

A.1 Proof of theorem 1

Clearly, Theorem 1 relies on the asymptotic linearity of the estimators $\hat{P}_{n_{p},hj}^{(p)}(0,t)$, $p=1,2$. This can be established by first showing the asymptotic linearity of the Nelson–Aalen estimators of the cumulative transition intensities and then by utilizing the functional delta method (Van der Vaart, 2000). The steps to achieve this utilize conditions C1 and C3–C5 and arguments similar to those used in the proof of Theorem 2 of Bakoyannis et al. (2019). After this analysis it can be shown that

[TABLE]

with

[TABLE]

By Lemma 1 in the supplementary material of Bakoyannis et al. (2019) and arguments similar to those use in the proof of Theorem 2 of the same source, it follows that the influence functions $\gamma_{ihj}^{(p)}(0,t)$, $p=1,2$, belong to $P$-Donsker classes of functions. Now, it is not hard to see that under the null hypothesis and by conditions C2 and C6

[TABLE]

Finally, the statement of Theorem 1 follows as a result of the usual central limit theorem and the independence between the two terms, as a consequence of the fact that the two samples are independent.

A.2 Proof of theorem 2

Due to the asymptotic linearity of the transition probability estimators $\hat{P}_{n_{p},hj}^{(p)}(0,t)$, for $p=1,2$, as argued in the proof of Theorem 1, along with conditions C2 and C6, it follows that

[TABLE]

Now, by the Donsker property of the class of functions $\{\gamma_{hj}^{(p)}(0,t):p=1,2,t\in[0,\tau]\}$ and the uniform boundedness of the class of fixed functions $\{W_{hj}(t):t\in[0,\tau]\}$, it follows that $\{W_{hj}(t)\gamma_{hj}^{(p)}(0,t):p=1,2,t\in[0,\tau]\}$ is also a $P$-Donsker class. Therefore, by the independence between the two samples, it follows that

[TABLE]

where $\mathbb{G}_{1hj}$ and $\mathbb{G}_{2hj}$ are two independent tight zero-mean Gaussian processes with covariance functions

[TABLE]

Now, define

[TABLE]

where $\xi^{(p)}$, $p=1,2$, are independent random draws from $N(0,1)$. By the Donsker property of the class $\{W_{hj}(t)\gamma_{hj}^{(p)}(0,t):p=1,2,t\in[0,\tau]\}$, for $h,j\in\mathcal{I}$, and the conditional multiplier central limit theorem (Van Der Vaart and Wellner, 1996) it follows that

[TABLE]

conditionally on the observed data. Therefore

[TABLE]

conditionally on the observed data. Now it remains to argue that $\sup_{t\in[0,\tau]}|\hat{B}_{hj}(t)-\bar{B}_{hj}(t)|\equiv\|\hat{B}_{hj}(t)-\bar{B}_{hj}(t)\|_{\infty}=o_{p}(1)$, unconditionally on the observed data. By the triangle inequality it follows that

[TABLE]

By similar calculations to those in the proof of Theorem 2 in Bakoyannis et al. (2019) and conditions C1-C6 it follows that both normed terms in right side the above inequality are $o_{p}(1)$. This concludes the proof of Theorem 2.