A Reproducing Kernel Hilbert Space log-rank test for the two-sample problem

TL;DR

This paper introduces a novel two-sample test for right-censored data using a supremum of log-rank tests indexed by a reproducing kernel Hilbert space, offering robustness and simplicity.

Contribution

It proposes a new RKHS-based log-rank test that is exact, easy to evaluate, and connects with existing methods, enhancing robustness against various alternatives.

Findings

The test is omnibus for certain RKHS families.

The method provides an exact evaluation of the test statistic.

Empirical evaluation demonstrates its effectiveness.

Abstract

Weighted log-rank tests are arguably the most widely used tests by practitioners for the two-sample problem in the context of right-censored data. Many approaches have been considered to make weighted log-rank tests more robust against a broader family of alternatives, among them, considering linear combinations of weighted log-rank tests, and taking the maximum among a finite collection of them. In this paper, we propose as test statistic the supremum of a collection of (potentially infinite) weight-indexed log-rank tests where the index space is the unit ball in a reproducing kernel Hilbert space (RKHS). By using some desirable properties of RKHSs we provide an exact and simple evaluation of the test statistic and establish connections with previous tests in the literature. Additionally, we show that for a special family of RKHSs, the proposed test is omnibus. We finalise by performing an empirical evaluation of the proposed methodology and show an application to a real data scenario. Our theoretical results are proved using techniques for double integrals with respect to martingales that may be of independent interest.

Tables

Table 1. Table 1: Significance Level at α = 5 % 𝛼 percent 5 \alpha=5\% . n c subscript 𝑛 𝑐 n_{c} stands for the number of data points of group c 𝑐 c , and cen c subscript cen 𝑐 \text{cen}_{c} indicates the percentage of censored data.

				Log-Rank		Projection		Pearson-type		SEK
$n_0$	$n_1$	${\text{cen}}_0$	${\text{cen}}_1$	LRP	LRC	P2W	P4W	Per4	Per5	SEK
30	30	10%	10%	4.5	5	5.1	4.7	6.1	5.7	5
		10%	30%	5.4	4.4	4.5	4.1	5	5.1	4.2
		30%	30%	4.3	4.7	4.2	4.4	4.7	4.4	4.2
30	100	10%	10%	3.9	5.7	5	5.1	5.2	4.9	5.2
		10%	30%	6.7	4.6	7.1	5.1	5.2	4.5	4.7
		30%	10%	5	7.3	6.1	4.9	5.9	5.6	5
		30%	30%	3.9	5.1	3.8	3.9	4.4	4.8	3.8
100	100	10%	10%	4.7	4.2	5.2	5.9	4.9	4.7	4.3
		10%	30%	5	5.9	5.4	6.6	6.1	5.8	6.1
		30%	30%	5.1	5.3	5.5	4.7	5.1	4.8	4.8

Table 2. Table 2: p 𝑝 p -values obtained (using the Wild Bootstrap approach) for the comparison of the treatments ‘Chemotherapy’ and ‘Chemotherapy+Radiation’ for the GTSG data.

	SEK	Per4	Per5	LRP	LRC	P2W	P4W
p-values	0.0053	0.0151	0.0222	0.2531	0.0011	0.0051	0.0228

Equations

$$\left\{\text{mathds}{1}_{\{X_i\le s,{\Delta }_i=0\}},\text{mathds}{1}_{\{X_i\le s,{\Delta }_i=1\}}:0\le s\le x,i\in [n]\right\},$$

$${\mathrm{LR}}_n(\omega )=\frac{n}{n_0{n}_1}{\int }_0^{{\tau }_n}\omega (\hat{F}(x-))L(x)\left(d{\hat{\Lambda }}_0(x)-d{\hat{\Lambda }}_1(x)\right),$$

$$\sqrt{\frac{n_0{n}_1}{n}}{\mathrm{LR}}_n(\omega )$$

$${\sigma }^2$$

$$\Lambda (x;\theta ,\omega )$$

$$U(0)$$

$$\hat{U}(0)$$

$$Z_n$$

$$K(x,y)=⟨K_x,K_y{⟩}_{\mathcal{H}}.$$

$$\mu \to {\int }_{[0,1)}K(y,\cdot )\mu (dy)\in \mathcal{H},$$

$$Z_n=\underset{\omega \in \mathcal{H},\Vert \omega {\Vert }_{\mathcal{H}}^2\le 1}{\sup }{\mathrm{LR}}_n(\omega {)}^2,$$

$$Z_n$$

$$Z_n={\left(\frac{n}{n_0{n}_1}\right)}^2{\mathbf{V}}^{⊺}\hat{\mathbf{K}}\mathbf{V},$$

$${\varphi }_{\nu }(\cdot )$$

$${\varphi }_0^n(\cdot )=\frac{n}{n_0{n}_1}{\int }_0^{{\tau }_n}K(\hat{F}(y-),\cdot )L(y)d{\hat{\Lambda }}_0(y),\text{and}{\varphi }_1^n(\cdot )=\frac{n}{n_0{n}_1}{\int }_0^{{\tau }_n}K(\hat{F}(y-),\cdot )L(y)d{\hat{\Lambda }}_1(y),$$

$${\nu }_0^n((s,t])=\frac{n}{n_0{n}_1}{\int }_s^{t}L(x)d{\hat{\Lambda }}_0(x),$$

$${\mathrm{LR}}_n(\omega )$$

$$Z_n=\Vert {\varphi }_0^n-{\varphi }_1^n{\Vert }_{\mathcal{H}}^2.$$

$${\mathrm{LR}}_n(\omega )$$

$$=\frac{n}{n_0{n}_1}{\int }_0^{{\tau }_n}⟨w,K(\hat{F}(x-),\cdot ){⟩}_{\mathcal{H}}L(x)\left(d{\hat{\Lambda }}_0(x)-d{\hat{\Lambda }}_1(x)\right)$$

$$={⟨w,\frac{n}{n_0{n}_1}{\int }_0^{{\tau }_n}K(\hat{F}(x-),\cdot )L(x)\left(d{\hat{\Lambda }}_0(x)-d{\hat{\Lambda }}_1(x)\right)⟩}_{\mathcal{H}}$$

$$={⟨w,{\varphi }_0^n-{\varphi }_1^n⟩}_{\mathcal{H}}.$$

$$Z_n$$

$$\Vert {\varphi }_0^n-{\varphi }_1^n{\Vert }_{\mathcal{H}}^2$$

$$={⟨\frac{n}{n_0{n}_1}{\int }_0^{{\tau }_n}K(\hat{F}(x-),\cdot )L(x)\left(d{\hat{\Lambda }}_0(x)-d{\hat{\Lambda }}_1(x)\right),\frac{n}{n_0{n}_1}{\int }_0^{{\tau }_n}K(\hat{F}(y-),\cdot )L(y)\left(d{\hat{\Lambda }}_0(y)-d{\hat{\Lambda }}_1(y)\right)⟩}_{\mathcal{H}}$$

$$={\left(\frac{n}{n_0{n}_1}\right)}^2{\int }_0^{{\tau }_n}{\int }_0^{{\tau }_n}{⟨K(\hat{F}(x-),\cdot ),K(\hat{F}(y-),\cdot )⟩}_{\mathcal{H}}L(x)L(y)\left(d{\hat{\Lambda }}_0(x)-d{\hat{\Lambda }}_1(x)\right)\left(d{\hat{\Lambda }}_0(y)-d{\hat{\Lambda }}_1(y)\right)$$

$$={\left(\frac{n}{n_0{n}_1}\right)}^2{\int }_0^{{\tau }_n}{\int }_0^{{\tau }_n}K(\hat{F}(x-),\hat{F}(y-))L(x)L(y)\left(d{\hat{\Lambda }}_0(x)-d{\hat{\Lambda }}_1(x)\right)\left(d{\hat{\Lambda }}_0(y)-d{\hat{\Lambda }}_1(y)\right).$$

$$Z_n={\mathrm{LR}}_n(\omega {)}^2.$$

$$Z_n=\sum _{j=0}^{k-1}{\mathrm{LR}}_n{\left(\omega \text{mathds}{1}_{\{I_j\}}\right)}^2=\sum _{j=0}^{k-1}{\left({\int }_{I_j}\omega (\hat{F}(x-))L(x)\left(d{\hat{\Lambda }}_0(x)-d{\hat{\Lambda }}_1(x)\right)\right)}^2.$$

$$Z_n$$

$$K(x,y)=\sum _{j=0}^{k-1}\frac{\omega (x)\omega (y)}{{\sigma }_j(\omega {)}^2}\text{mathds}{1}_{I_j\times I_j}(x,y),$$

$$\nu (dx)$$

Full text

Key words: Survival Analysis, Right-Censored Data, Reproducing Kernel Hilbert Space, Log-rank Test, two-sample tests.

1 Introduction

Two-sample testing is a classical problem in the context of survival data. For instance, in a clinical trial, two-sample tests can be used to compare different treatments when the survival times of patients are censored. Within the context of right-censored data, the classical log-rank test, first introduced by Mantel (1966) and Peto and Peto (1972), is the most widely used test among practitioners. A well-known property of the classical log-rank test is that it is the most powerful test under the assumption of proportional-hazard alternatives. This result can be deduced by noticing that the log-rank test statistic coincides with the score test statistic when the true cumulative hazard function belongs to the model $\left\{\Lambda_{\theta},\theta\in\Theta\right\}$, where $\Lambda_{\theta}(t)=\int_{0}^{t}e^{\theta}d\Lambda_{0}(x)$, $\Lambda_{0}$ denotes the cumulative hazard function under the null hypothesis, and $\theta\in\Theta\subseteq\mathds{R}$ is the parameter of the model. While the classical log-rank test is optimal for proportional-hazard alternatives, it can have a substandard behaviour when the true cumulative hazard function cannot be expressed in terms of $\{\Lambda_{\theta},\theta\in\Theta\}$.

In order to broaden the power of the classical log-rank test to other families of alternatives, researchers have introduced and studied different variants of weighted log-rank tests (Tarone and Ware, 1977; Gill, 1980; Harrington and Fleming, 1982; Bagdonavicius et al., 2010; Andersen et al., 2012). We refer the reader to Chapter 7 of Klein and Moeschberger (2006) for a general discussion and comparison of different weighted log-rank tests approaches. In the simplest case, each weighted log-rank test with weight function $\omega$ can be written as a score test for alternatives in the model $\left\{\Lambda_{\theta,\omega}(t),\theta\in\Theta\right\}$, where $\Lambda_{\theta,\omega}(t)=\int_{0}^{t}e^{\theta\omega(F_{0}(x))}d\Lambda_{0}(x)$ and $\omega$ is fixed. Then, similarly to the above, it can be deduced that each weighted log-rank test is the most powerful test for the null hypothesis $H_{0}:\Lambda_{1}=\Lambda_{0}$ under the assumption that $\Lambda_{1}$ can be expressed in terms of $\left\{\Lambda_{\theta,\omega},\theta\in\Theta\right\}$, for some $\theta\neq 0$. If the true cumulative hazard function cannot be expressed in terms of this parametric model, there are no guarantees at all for the behaviour of the weighted log-rank test. Indeed, it may happen that we observe pathological cases in which the test has zero asymptotic power for specific alternatives. This is the case of the classical log-rank test, which is recovered by choosing $\omega=1$, in the setting of crossing-hazard alternatives.

In an attempt to overcome the previously described limitations of weighted log-rank tests, researchers have considered two natural approaches to improve their performance: the first approach considered is to learn the weight function from the data, which defines a weighted log-rank test with an adaptive weight function, and the second approach consists on combining several log-rank tests into a single test statistic. Both approaches avoid making the strong parametric assumption that a particular model is true. The first approach has been discussed by Lai et al. (1991), Yang et al. (2005) and Yang and Prentice (2010), and the second approach, which is the focus of this paper, has been studied by several authors. Particularly, Tarone (1981), Kössler (2010), and Garès et al. (2015) combined weighted log-rank tests by considering the maximum of a finite collection of them, and Kosorok and Lin (1999) proposed the supremum of an infinite collection of weight-indexed log-rank tests, with weights belonging to a particular space of functions. The latter approach is very close to what we propose in this paper, however, the test statistic proposed by Kosorok and Lin (1999) lacks of an analytically tractable expression, which forces the authors to rely on a Monte-Carlo approximation of it, whereas our approach, by an appropriate selection of the space of functions, obtains a simple expression for our test statistic, leading to a simple testing procedure.

While the focus of this paper are weighted log-rank test and the combination of them, there are other test-approaches for the two-sample problems that are worth mentioning. One example of these approaches are the weighted Kaplan-Meier estimators (Pepe and Fleming, 1989), which are parametrised on a weight function, and they can also be combined into one single statistic (Shen and Cai, 2001). A different approach is the so-called empirical likelihood approach (Zhou, 2016), which leads to very interesting test statistics that can also be combined into a single statistic (Bathke et al., 2009). Somewhat classical options are the Cramer von-Mises and the Kolmovorov-Smirnov test statistics, that have been deeply studied in many settings, including Survival Analysis (Koziol and Green, 1976; Koziol and Yuh, 1982).

In this paper, we consider the supremum of a potentially infinite collection of weighted log-rank tests with weights belonging to the unit ball of a reproducing kernel Hilbert space of functions. Then, we propose to test the null hypothesis $H_{0}:\Lambda_{1}=\Lambda_{0}$, by checking if $\sup_{\omega\in\mathcal{H}:\|\omega\|_{\mathcal{H}}\leq 1}\text{LR}_{n}(\omega)\approx 0$, where $\text{LR}_{n}(\omega)$ denotes a weighted log-rank test with weight function $\omega$, and $\mathcal{H}$ is a reproducing kernel Hilbert space of functions. Clearly, if $\Lambda_{1}$ can be expressed in terms of $\{\Lambda_{\theta,\omega},\theta\in\Theta,\omega\in B_{1}(\mathcal{H})\}$ for some $\theta\neq 0$ and some $\omega\in B_{1}(\mathcal{H})$ (where $B_{1}(\mathcal{H})$ denotes the unit ball of $\mathcal{H}$), our test statistic should be strictly greater than zero, in which case we decide against the null hypothesis. Also, notice that, the larger the space of functions $\mathcal{H}$, the more models we test at the same time.

Choosing weights $\omega$ in the unit ball of a reproducing kernel Hilbert space $\mathcal{H}$ is the key step in our testing procedure, as we will show that, by doing so, our test statistic $\sup_{\omega\in\mathcal{H}:\|\omega\|_{\mathcal{H}}\leq 1}\text{LR}_{n}(\omega)$ can be evaluated exactly by a straightforward computation. This result follows from the good properties of reproducing kernel Hilbert spaces. Particularly, we will use two properties of reproducing kernel Hilbert spaces: i) they are uniquely characterised by a kernel function $K$, i.e., a symmetric and positive-definite function, and ii) they satisfy the so-called reproducing kernel property, which states that $\omega(x)=\langle\omega,K(x,\cdot)\rangle_{\mathcal{H}}$ holds for any function $\omega\in\mathcal{H}$ and any point $x$. By using this property we will show that $\sup_{\omega\in\mathcal{H}:\|\omega\|_{\mathcal{H}}}\text{LR}_{n}(\omega)^{2}=\sum_{i=1}^{n}\sum_{j=1}^{n}K(x_{i},x_{j})c_{i}c_{j}$, for some points $x_{i}$ and constants $c_{j}$ (depending on the data). This type of kernel estimators have been previously addressed by Neuhaus (1987) in connection with infinite sums of linear rank estimators (such as the log-rank) in the uncensored setting. Since, in practice, our test statistic is a quadratic form based on the kernel $K$, our method can be described just in terms of $K$ avoiding the need to describe a Hilbert space. Describing kernel functions with interesting properties and useful interpretation has been a topic of extensive research, especially in the Machine Learning community. A compendium of popular kernels used in applications can be found in Sousa (2010).

While we can easily establish connections between our method and approaches based on weighted log-rank tests, we will show that an alternative interpretation of our test statistic allows us to connect our approach with the kernel mean embedding testing approach. In particular, we will show that our test statistic $\sup_{\omega\in\mathcal{H}:\|\omega\|_{\mathcal{H}}\leq 1}\text{LR}_{n}(\omega)$ can be alternatively written as the norm of the difference of particular embeddings of our two samples into $\mathcal{H}$. The idea of comparing probability distributions/datasets by embedding them into a reproducing kernel Hilbert space of functions has been extensively studied in the uncensored setting by researchers in Statistics and Machine Learning (Berlinet and Thomas-Agnan, 2004; Smola et al., 2007), however, it seems this idea has yet to percolate into the Survival Analysis community and, up to the best of our knowledge, it has only been considered by Fernandez and Gretton (2019) for the goodness-of-fit problem under right-censored data. In the uncensored case, two-sample tests were proposed and studied by Gretton et al. (2012). In this work, the authors embed two different samples into a reproducing kernel Hilbert space by considering the so-called “kernel mean embedding” of the empirical distribution of each sample. Then, they compute the difference between these two samples by computing the distance (induced by the norm) of the embeddings into $\mathcal{H}$. This idea can be straightforwardly applied to right-censored data by considering the kernel mean embedding of the Kaplan-Meier estimator instead of the empirical distribution, which leads to test statistics that are Kaplan-Meier $V$-statistics. These type of statistics were studied by Bose and Sen (2002) and by Fernández and Rivera (2020), and include as particular cases well-known statistics such as the Cramer von-Mises test statistic. As it was pointed out by Fernández and Rivera (2020) such an approach is not suitable for censored data as Kaplan-Meier $V$-statistics are only reliable when the amount of censored observations is rather small when compared to the total amount of data, and even in that case, the test statistic may not be data-efficient.

In this work, we study several asymptotic properties of our test statistic under both, the null and alternative hypothesis. Under the null hypothesis, we find the limiting distribution of our test statistic by approximating it by a degenerate V-statistic. Under the alternative hypothesis, we prove that under reasonable conditions, our test has asymptotic power tending to one. While the asymptotic null distribution is known, in most cases, its parameters are hard to compute. Thus, in order to implement our test, we propose a Wild Bootstrap approximation of the null distribution and prove that the Wild Bootstrap statistic converges in distribution to our test statistic under the null hypothesis, giving us a correct Type-I error. Finally, we show that particular instances of our test statistic recover some existing testing approaches studied in the literature, which suggests our testing procedure is a natural generalisation of them. Examples of tests which our testing procedure recovers are: weighted log-rank tests, Pearson-type tests (Akritas, 1988) and projection-type tests (Brendel et al., 2014; Ditzhaus and Friedrich, 2018).

Besides theoretical guarantees, we provide an extensive empirical evaluation of our method in two important data scenarios: proportional hazard functions and time-dependent hazard functions, including Weibull and periodic hazard alternatives. In our experiments we demonstrate that our method has a good performance in a wide range of problems, which include sample sizes from small to large, and different censoring percentages. Our experiments show that finite-dimensional reproducing kernel Hilbert spaces tend to have an overall better performance in problems with fewer observations or a simpler hazard function. While more complex reproducing kernel Hilbert spaces (i.e., infinite-dimensional spaces) still show a good performance in these simpler problems, our experiments show that they are better suited for larger datasets and more complex hazard functions. We also provide a real-data evaluation of our method, trying different kernel functions, and we compare the results with those obtained by Ditzhaus and Friedrich (2018).

The structure of the paper is as follows. In Section 2, we introduce some standard notation used in Survival Analysis and describe weighted log-rank tests. In Section 3 we introduce the essential background knowledge about reproducing kernel Hilbert spaces (RKHS) and formally introduce our test statistic. In Section 4 we study asymptotic properties of our test. Later, in Section 5, we proceed to explain how to implement a Wild Bootstrap approach. Sections 6 and 7 are devoted to empirical studies using simulated data and real data, respectively.

2 Survival analysis background

2.1 General notation

We establish some general notation that will be used throughout the paper. We denote $\mathds{R}_{+}=(0,\infty)$. Let $f:\mathds{R}_{+}\to\mathds{R}_{+}$ be an arbitrary right-continuous function, then we define $f(x-)=\lim_{h\to 0}f(x-|h|)$. In this work we make use of standard asymptotic notation (Janson, 2011), e.g., $O_{p}$, $o_{p}$, $\Theta_{p}$, etc., with respect to the number of sample points $n$. In order to avoid large parentheses, we write $X=O_{p}(1)Y$ instead of $X=O_{p}(Y)$, especially if the expression for $Y$ is very long. Given a sequence of stochastic processes $(W_{n}(x):x\in\mathcal{X})$, depending on the number of observations $n$, and a function $f(x)\geq 0$, we say that $W_{n}(x)=O_{p}(1)f(x)$ uniformly for all $x\in A_{n}$, if and only if $\sup_{x\in A_{n}}|W_{n}(x)|/f(x)=O_{p}(1)$, where $A_{n}\subseteq\mathcal{X}$ is a set that may depend on $n$, and we use the convention $0/0=0$. Lastly, in this work the integral symbol $\int_{a}^{b}$ means integration over $(a,b]$ unless we state otherwise.

2.2 Right-censored data

Our data corresponds to right-censored observations belonging to two groups/classes, namely group [math] and group 1. We assume that the total number of observations is $n$, and that $n=n_{0}+n_{1}$, where $n_{c}$ denotes the number of observations in group $c\in\{0,1\}$. We use $[n]$ to denote the set $\{1,\ldots,n\}$. For the sake of asymptotic analysis, we assume that we have no vanishing groups, that is, $n_{0}/n\to\eta$ and $n_{1}/n\to(1-\eta)$ with $\eta\in(0,1)$ as $n$ tends to infinity.

Right-censored datasets are commonly observed in triples $((X_{i},\Delta_{i},c_{i}))_{i\in[n]}$, where $X_{i}=\min\{T_{i},C_{i}\}$ is an observed time, defined as the minimum between a survival time of interest $T_{i}$ and a censoring time $C_{i}$, the variable $\Delta_{i}=\mathds{1}_{\{X_{i}=T_{i}\}}$ is an indicator of whether we actually observe the survival time of interest, and $c_{i}\in\{0,1\}$ is an associated covariate that denotes the group membership of the $i$-th observation. In this work we assume that all triples are mutually independent and that the triples have the same distribution within each group. Additionally, we assume non-informative censoring, meaning that, given the covariate $c_{i}$, the survival time $T_{i}$ is independent of the censoring time $C_{i}$.

We denote by $F_{c}$, $G_{c}$, and $H_{c}$, the respective conditional distribution functions of the random variables $T_{i}$, $C_{i}$, and $X_{i}$, given that the covariate $c_{i}$ takes the value $c\in\{0,1\}$. Notice that $1-H_{c}=(1-F_{c})(1-G_{c})$ due to the non-informative censoring assumption. For simplicity of exposition, we assume that the distribution functions $F_{c}$ and $G_{c}$ are continuous distributions on $\mathds{R}_{+}$, however, all the methods of this paper can be extended to general probability distributions on $\mathds{R}$, as our results are based on counting processes arguments that have been developed in full generality. We denote by $S_{c}(x)=1-F_{c}(x)$ and by $\Lambda_{c}(x)=\int_{(0,x]}S_{c}(s)^{-1}dF_{c}(s)$, the so-called survival and cumulative hazard functions associated with $F_{c}$, respectively, and we denote by $\widehat{S}_{c}$ the Kaplan-Meier estimator of $S_{c}$, and by $\widehat{\Lambda}_{c}$ the Nelson-Aalen estimator of $\Lambda_{c}$. Notice that for the computation of $\widehat{S}_{c}$ and $\widehat{\Lambda}_{c}$ we only use triples $(X_{i},\Delta_{i},c_{i})$ satisfying $c_{i}=c$.

While in principle we assume that the covariates $c_{i}$ are deterministic, it is also possible to consider them as independent and identically distributed (i.i.d.) samples from a Bernoulli distribution with mean $\eta$ (notice that in this case $n_{0}$ is a random variable satisfying $n_{0}/n\to\eta$ almost surely). These two possible assumptions originate two models: the deterministic covariates model and the random covariates model. For our analysis, it will be convenient to use the random covariates model as, under this assumption, the triples $((X_{i},\Delta_{i},c_{i}))_{i\in[n]}$ are i.i.d. which simplifies the asymptotic analysis, nevertheless, as we will show in Appendix B.3, both models are asymptotically equivalent. Also, in practice, there is no difference between the two models when implementing our testing procedure.

In this paper, we use the notion of pooled data and pooled distributions. The pooled data corresponds to the original data after ignoring/removing the covariates $(c_{i})_{i\in[n]}$, that is, the pooled data is $((X_{i},\Delta_{i}))_{i\in[n]}$. We use the term pooled distributions to refer to the distributions of a randomly selected data-point. Particularly, under the random covariates model, the pooled distributions associated with the observed times $(X_{i})_{i\in[n]}$ and the survival times $(T_{i})_{i\in[n]}$, are the marginal distributions of $X_{i}$ and $T_{i}$, respectively, given by $H=\eta H_{0}+(1-\eta)H_{1}$ and $F=\eta F_{0}+(1-\eta)F_{1}$. Under the deterministic covariates model, the pooled distributions associated with the observed times and the survival times correspond to $(n_{0}/n)H_{0}+(n_{1}/n)H_{1}$ and $(n_{0}/n)F_{0}+(n_{1}/n)F_{1}$, respectively. Notice that, since $(n_{0}/n)H_{0}+(n_{1}/n)H_{1}\to H$ and $(n_{0}/n)F_{0}+(n_{1}/n)F_{1}\to F$ as $n$ grows to infinity, the pooled distributions of both models are asymptotically equivalent. Additionally, we denote by $S(x)=1-F(x)$ and by $\Lambda(x)=\int_{(0,x]}S(s)^{-1}dF(s)$ the survival and cumulative hazard functions associated with the pooled distribution $F$, and we use the Kaplan-Meier estimator $\widehat{S}$ and the Nelson-Aalen estimator $\widehat{\Lambda}$ to approximate $S$ and $\Lambda$, respectively. Notice that $\widehat{S}$ and $\widehat{\Lambda}$ are computed using the pooled-data.

2.3 Counting processes

In this work, we use the standard counting processes notation used in Survival Analysis. We define the individual, class and pooled counting processes by $N^{i}(x)=\Delta_{i}\mathds{1}_{\{X_{i}\leq x\}}$ for $i\in[n]$, $N_{c}(x)=\sum_{i=1}^{n}\mathds{1}_{\{c_{i}=c\}}N^{i}(x)$ for $c\in\{0,1\}$, and $N(x)=\sum_{i=1}^{n}N^{i}(x)=N_{0}(x)+N_{1}(x)$, respectively. Similarly, we define the individual, class and pooled risk functions by $Y^{i}(x)=\mathds{1}_{\{X_{i}\geq x\}}$, $Y_{c}(x)=\sum_{i=1}^{n}\mathds{1}_{\{c_{i}=c\}}Y^{i}(x)$ and $Y(x)=\sum_{i=1}^{n}Y^{i}(x)=Y_{0}(x)+Y_{1}(x)$, respectively. By using the previous notation, we write the Nelson-Aalen estimator of $d\Lambda_{c}$ as $d\widehat{\Lambda}_{c}(x)=dN_{c}(x)/Y_{c}(x)$, and the pooled Nelson-Aalen estimator as $d\widehat{\Lambda}(x)=dN(x)/Y(x)$. To avoid ambiguities we assume that $0/0=0$, e.g., $d\widehat{\Lambda}(x)=0$ if $Y(x)=0$. Finally, we define the process $L(x)=Y_{0}(x)Y_{1}(x)/Y(x)$, which appears frequently in our results, and we define $\tau_{n}=\max\{X_{i}:i\in[n]\}$ and $\tau=\sup\{x:1-H(x)>0\}$. Observe that $\tau_{n}\to\tau$ a.s. and that $N^{i}(x)=N^{i}(\tau_{n})$, $N_{c}(x)=N_{c}(\tau_{n})$ and $N(x)=N(\tau_{n})$ for any $x\geq\tau_{n}$. Also, notice that all our observed times, $X_{i}$, are less than $\tau$ as they are generated from continuous distributions.

We assume that all random variables are defined on a common filtrated probability space $(\Omega,\mathcal{F},(\mathcal{F}_{x})_{x\geq 0},\mathds{P})$, where the sigma-algebra $\mathcal{F}_{x}$ is generated by

[TABLE]

and the $\mathds{P}$-null sets of $\mathcal{F}$. Under the deterministic covariates model, we define the individual $(\mathcal{F}_{x})$-martingales associated to each data point, $(X_{i},\Delta_{i},c_{i})$, as $M^{i}(x)=N^{i}(x)-\int_{(0,x]}Y^{i}(s)d\Lambda_{c_{i}}(s)$. Similarly, under the deterministic covariates model, we define the class and pooled martingales by $M_{c}(x)=\sum_{i=1}^{n}\mathds{1}_{\{c_{i}=c\}}M^{i}(x)$ and $M(x)=\sum_{i=1}^{n}M^{i}(x)=M_{0}(x)+M_{1}(x)$, respectively, where $c\in\{0,1\}$. Notice that under the random covariates model, $M^{i}(x)$, $M_{c}(x)$ and $M(x)$ are $(\mathcal{F}_{x})$-martingales when conditioning on $(c_{i})_{i\in[n]}$ (alternatively, we can include the sigma algebra associated with $(c_{i})_{i\in[n]}$ to generate $\mathcal{F}_{x}$). We denote by $\langle M^{i}\rangle(x)$ and $[M^{i}](x)$, respectively, the predictable and quadratic variation processes associated with $M^{i}(x)$, which, in the context of continuous survival and censoring times, are given by $d\langle M^{i}\rangle(x)=Y^{i}(x)d\Lambda_{c_{i}}(x)$ and $d[M^{i}](x)=dN^{i}(x)$. For further information about counting processes and their applications in Survival Analysis we refer the reader to Fleming and Harrington (1991).

Notice that, since we consider continuous survival times, we can estimate $F(x)$ by either $\widehat{F}(x)$ or $\widehat{F}(x-)$. The main advantage of using the latter estimator is that $(\widehat{F}(x-))_{x\geq 0}$ is left-continuous, and thus, it is a predictable process with respect to the filtration $(\mathcal{F}_{x})_{x\geq 0}$.

2.4 The log-rank estimator

The weighted log-rank statistic is defined by

[TABLE]

where the function $\omega:[0,1)\to\mathds{R}$ is referred to as weight function. In the previous equation recall that $L(x)=Y_{0}(x)Y_{1}(x)/Y(x)$, $\widehat{F}$ is the pooled Kaplan-Meier estimator of the pooled distribution function $F$, and $\widehat{\Lambda}_{0}$ and $\widehat{\Lambda}_{1}$ are the class Nelson-Aalen estimators of the cumulative hazard functions $\Lambda_{0}$ and $\Lambda_{1}$, respectively.

Gill (1980) studied the asymptotic behaviour of the weighted log-rank statistic, $\operatorname{LR}_{n}(\omega)$, under the null hypothesis $H_{0}:F_{0}=F_{1}$, obtaining that

[TABLE]

where

[TABLE]

Notice that, even under the null hypothesis, the censoring distributions $G_{0}$ and $G_{1}$ are not necessarily equal, hence the expression given for $\sigma^{2}$ cannot be simplified. Also, in most cases, the censoring distributions are unknown, and thus it is not possible to evaluate $\sigma^{2}$ to find rejection regions. To carry-out a testing procedure, we can use the asymptotic null distribution given in Equation (2), replacing the asymptotic variance $\sigma^{2}$ by $\widehat{\sigma}^{2}=n/(n_{0}n_{1})\int_{0}^{\tau_{n}}\omega(\widehat{F}(x-))^{2}L(x)d\widehat{\Lambda}(x)$, which is an unbiased estimator of $\sigma^{2}$ under the null hypothesis. We refer the reader to Gill (1980) for more details.

It is a well-known fact that weighted log-rank tests relate to score tests through the choice of a specific family of alternatives. Let

[TABLE]

be a parametric model (indexed by $\theta$) for the cumulative hazard function, where $\omega:[0,1)\to\mathds{R}$ is a fixed continuous function, and $\Theta$ is chosen as an open subset of $\mathds{R}$ containing [math]. Then, under the assumption that $\Lambda_{1}(x)=\Lambda(x;\theta_{1},\omega)$ for some $\theta_{1}\in\Theta$, testing the null hypothesis, $H_{0}:\Lambda_{0}=\Lambda_{1}$, is equivalent to test $H_{0}:\theta=0$.

We can test $H_{0}:\theta=0$ by using a score test. The score statistic is computed in terms of the score function defined by $U(\theta)=\partial/\partial\theta\log L_{n}(\theta;\omega)$, where $L_{n}(\theta;\omega)$ is the likelihood function under the model of Equation $\eqref{eqn:familyAternal}$. For the Goodness-of-Fit problem, where $F_{0}$ and $\Lambda_{0}$ are known, the score statistic is given by

[TABLE]

In the Two-Sample problem, $\Lambda_{0}$ and $F_{0}$ are unknown, but they can be estimated using the pooled Nelson-Aalen and Kaplan-Meier estimators, obtaining

[TABLE]

Then, a simple comparison shows that $n/(n_{0}n_{1})\widehat{U}(0)=\operatorname{LR}_{n}(\omega)$, from which we deduce the relationship between the weighted log-rank test defined in Equation (1) and the score test associated to the parametric model defined in Equation (3). By the Neyman-Pearson Lemma, we deduce that the weighted log-rank test is the most powerful test for small departures from the null, that is, when $\theta\to 0$. Unfortunately, if the model in Equation (3) is misspecified, little can be said about the performance of the test, indeed, it is well-known that, in some cases, the weighted log-rank statistic yields a test with asymptotically zero-power.

3 An RKHS approach to the log-rank test

A standard approach to broaden the power of weighted log-rank tests, in order to achieve a more robust behaviour across a larger class of alternatives, is to combine several weighted log-rank tests into a single statistic. In this way, if one of the weight functions completely fails to differentiate between the null and alternative hypotheses, we can still rely on the remaining weight functions to help us to discriminate, and thus increase the overall power of our testing procedure. Two interesting questions that arise from this approach are: How do we combine these weighted log-rank tests? Which and how many weight functions do we need to choose?

Selecting and combining weight functions efficiently for the problem at hand is very difficult as, in most cases, it requires that the user analyses the data in advance to select appropriate weight functions, e.g., it is usual to check if the hazard functions cross or/and if they show late/early departures. Searching for important features in the data may be very time-consuming and, moreover, it is always possible that there are features that do not translate into an appropriate weight function. Also, it might happen that the user overlooks some important relations.

Instead of relying on the expertise of a user to identify important features, we can consider as many weight functions as possible to account for the heterogeneity in the data. We can take this approach to the extreme by considering a potentially infinite family of weights. While this approach solves the problem of choosing weight functions, naively choosing a particular family of weights will lead to tests that, i) are hard to calibrate as we will need several data points to reach the correct Type-I error, ii) do not have an analytically tractable closed form for the test statistic, leading to tests that iii) are computationally expensive.

In order to overcome the previous difficulties, we consider weight functions belonging to the unit ball of a reproducing kernel Hilbert space $\mathcal{H}$, and propose as test statistic $Z_{n}$, which is defined by

[TABLE]

While a priori there is not a good reason to choose this particular family of weights, we will show that it has very nice features that translate into desirable properties of our testing procedure; among others, the test statistic has a simple closed form, allowing simple computations as well as an economic Wild Bootstrap implementation.

3.1 Reproducing kernel Hilbert spaces

A reproducing kernel Hilbert space (RKHS) is a Hilbert space of functions $\omega:[0,1)\to\mathds{R}$ satisfying that the evaluation $\omega\to\omega(x)$ is continuous for every fixed $x\in[0,1)$ (it is worth mentioning that we can replace $[0,1)$ by any space). By the Riesz representation theorem, for any $x\in[0,1)$, there exists a unique element $K_{x}\in\mathcal{H}$ such that $\omega(x)=\langle K_{x},\omega\rangle_{\mathcal{H}}$ for all $\omega\in\mathcal{H}$, which is known as the reproducing kernel property. Since $K_{x}\in\mathcal{H}$ for all $x\in[0,1)$, $K_{x}(y)=K_{y}(x)=\langle K_{x},K_{y}\rangle_{\mathcal{H}}$ holds for any $x,y\in[0,1)$. Then, this result allows us to define the so-called reproducing kernel $K:[0,1)^{2}\to\mathds{R}$ as

[TABLE]

From now on, in order to ease the notation, we write $K(x,\cdot)$ instead of $K_{x}(\cdot)$, even though the former induces a slight abuse of notation.

For every RKHS $\mathcal{H}$ with inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}}$ there exists a unique symmetric positive-definite reproducing kernel $K:[0,1)^{2}\to\mathds{R}$ satisfying Equation (5). Conversely, by the Moore-Aronszajn Theorem (Aronszajn, 1950), for any symmetric positive-definite kernel function $K:[0,1)^{2}\to\mathds{R}$, there exists a unique RKHS $\mathcal{H}$ for which $K$ is its reproducing kernel.

Given a finite signed measure $\mu$ on $[0,1)$, we define the kernel mean embedding of $\mu$ into $\mathcal{H}$ as

[TABLE]

where the previous integral has to be understood as a Bochner integral. A sufficient condition to guarantee the existence of such an embedding is that $\int_{[0,1)}\sqrt{K(x,x)}|\mu|(dx)<\infty$ (Gretton et al., 2012, Lemma 3). A kernel $K$ is said to be characteristic if the mean kernel embedding, defined in Equation (6), is injective on the space of probability distributions (i.e., distinct probability measures are embedded as different elements of $\mathcal{H}$). Furthermore, a continuous kernel $K$ is said to be $c$-universal if the mean kernel embedding is injective on the set of finite signed measures. Clearly, a continuous $c$-universal kernel is characteristic. Most of the standard kernels used in applications are $c$-universal, e.g., the Gaussian kernel, $K(x,y)=\exp\{-(x-y)^{2}/\sigma^{2}\}$, and the Ornstein–Uhlenbeck kernel, $K(x,y)=\exp\{-|x-y|/|\sigma|\}$. For more information about characteristic and $c$-universal kernels we refer the reader to Fukumizu et al. (2009), Muandet et al. (2017), Simon-Gabriel and Schölkopf (2018) and references therein.

3.2 An RKHS log-rank test

Recall that

[TABLE]

where $\mathcal{H}$ is an RKHS of functions $\omega:[0,1)\to\mathds{R}$ with reproducing kernel $K:[0,1)^{2}\to\mathds{R}$. In this section, we show that it is possible to find a closed-form expression for our test statistic, $Z_{n}$, in terms of the reproducing kernel $K$ associated with $\mathcal{H}$. This result is formally stated in the following theorem:

Theorem 1**.**

[TABLE]

By using the definition of the class Nelson-Aalen estimators, Equation (7) can be rewritten as

[TABLE]

which is a quadratic form with matrix-representation given by

[TABLE]

where $\mathbf{V}\in\mathds{R}^{n}$ is defined by $\mathbf{V}_{j}=L(X_{j})(-1)^{c_{j}}\Delta_{j}/Y_{c_{j}}(X_{j})$ and $\widehat{\mathbf{K}}\in\mathds{R}^{n}\times\mathds{R}^{n}$ is the matrix defined by $\widehat{\mathbf{K}}_{ij}=K(\widehat{F}(X_{i}),\widehat{F}(X_{j}))$.

In order to prove Theorem 1, we introduce an alternative embedding of the data, which is different to the kernel mean embedding of Equation (6). Given a measure $\nu$ (not necessarily a probability measure), we define the embedding $\phi$ of a finite signed measure $\nu$ into $\mathcal{H}$ by

[TABLE]

In practice, since the pooled distribution $F$ is unknown, we replace it by the pooled Kaplan-Meier estimator $\widehat{F}(x-)$ (recall that the survival time are continuous). Then, by using the previous definition, we introduce

[TABLE]

which are the corresponding embeddings of two empirical measures, $\nu^{n}_{0}$ and $\nu^{n}_{1}$, into $\mathcal{H}$, where $\nu^{n}_{0}$ and $\nu^{n}_{1}$ are defined for any $s<t$ by

[TABLE]

respectively. Notice that $\phi_{0}^{n}$ and $\phi_{1}^{n}$ are always well-defined, meaning that $\phi_{0}^{n}\in\mathcal{H}$ and $\phi_{1}^{n}\in\mathcal{H}$, since they are just finite sums of elements of $\mathcal{H}$ (observe that $K(\widehat{F}(y-),\cdot)\in\mathcal{H}$ for any fixed $y$).

By using the previous embeddings, we obtain an inner product representation of the weighted log-rank statistic, which will be used in the proof of Theorem 1.

Lemma 2** (Log-rank representation).**

For any $\omega\in\mathcal{H}$,

[TABLE]

and then

[TABLE]

Proof.

Since $\widehat{F}(x-)\in[0,1)$ for any $x\leq\tau_{n}$, the reproducing property yields $\omega(\widehat{F}(x-))=\langle\omega(\cdot),K(\widehat{F}(x-),\cdot)\rangle_{\mathcal{H}}$. Then, by the linearity of the inner product and integration,

[TABLE]

Finally, by taking supremum over the unit ball, we have that

[TABLE]

where the last equality holds since $\mathcal{H}$ is a Hilbert space. ∎

Proof of Theorem 1.

By Equation (10), we have $(\phi_{0}^{n}-\phi_{1}^{n})(\cdot)=\frac{n}{n_{0}n_{1}}\int_{0}^{\tau_{n}}K(\widehat{F}(x-),\cdot)L(x)(d\widehat{\Lambda}_{0}(x)-d\widehat{\Lambda}_{1}(x))$. Then, by using the linearity of the inner product and integration, and the reproducing kernel property, we deduce

[TABLE]

∎

3.3 Recovering existing tests

Our testing approach is based on fixing an RKHS of functions $\mathcal{H}$, which is done by fixing a kernel $K$. We show that for specific choices of the reproducing kernel $K$, we recover some previously known tests.

3.3.1 Weighted log-rank tests

In order to recover the classical weighted log-rank test $\operatorname{LR}_{n}(\omega)$, we set $K(x,y)=\omega(x)\omega(y)$. Notice that $K$ is symmetric and non-negative definite for any $\omega$. Then, by evaluating $Z_{n}$ using Equation (7), we obtain

[TABLE]

3.3.2 Pearson-type tests

Consider the partition of the interval $(0,1]$ given by $I_{j}=(j/k,(j+1)/k]$, where $j\in\{0,\ldots,k-1\}$ and $k>0$ is an integer. We recover Pearson-type testing approaches by choosing the kernel $K(x,y)=\sum_{j=0}^{k-1}\omega(x)\omega(y)\mathds{1}_{I_{j}\times I_{j}}(x,y),$ where $\omega:[0,1)\to\mathds{R}$ is a given weight function. Then, by evaluating $Z_{n}$ using Equation (7), we get

[TABLE]

Particularly, if we choose $\omega\equiv 1$, we recover

[TABLE]

which compares observed frequencies on $I_{j}$ between the two groups. We can go further and consider the normalised version of Pearson-type tests by choosing the kernel as

[TABLE]

where $\sigma_{j}(\omega)^{2}=\frac{n}{n_{0}n_{1}}\int_{I_{j}}\omega(\widehat{F}(x-))^{2}L(x)d\widehat{\Lambda}(x)$. Notice that in this case $K(x,y)$ is a random kernel, but it converges to a deterministic one since $\sigma_{j}(\omega)$ converges almost surely to a constant when $n$ tends to infinity.

3.3.3 Projection-type test

Projection-type tests were introduced by Brendel et al. (2014), and recently revisited by Ditzhaus and Friedrich (2018). Consider a finite number of weight functions $w_{1},\ldots,w_{k}$ such that $w_{i}\circ F\in L^{2}(\nu)$ for all $i\in\{1,\ldots,k\}$, where $\nu$ is the measure given by

[TABLE]

The projection-type testing approach considers the following test statistic

[TABLE]

where $\{\tilde{w}_{i}\}_{i=1}^{k}$ is an orthonormal base of the subspace $U$ of $L^{2}(\nu)$ generated by $w_{1}\circ F,\ldots,w_{k}\circ F$. We refer the reader to Brendel et al. (2014) for a detailed explanation. Brendel et al. (2014) and Ditzhaus and Friedrich (2018) recommend using weights with some meaning, for instance, $\omega(x)=x-1/2$ is used to detect a cross between two hazard functions around the median of the pooled survival time distribution, and $\omega(x)=x^{p}(1-x)^{q}$ is used to detect early and/or late differences between the hazard functions, depending on the parameters $p$ and $q$. Nevertheless, observe that, in terms of projections, the meaning of the functions does not matter as, for instance, projecting over the subspace generated by $\{1,x-1/2,x(1-x),x^{2}(1-x)\}$ is the same as projecting over the subspace generated by $\{1,x,x^{2},x^{3}\}$.

Our approach can be seen as a natural generalisation of the previous method. Indeed, we can recover the previous test statistic by choosing the kernel

[TABLE]

To compute the orthonormal basis $\{\tilde{\omega}_{i}\}_{i=1}^{k}$, we just need to compute the $k\times k$ matrix $\boldsymbol{P}$, where $\boldsymbol{P}_{ij}=\langle w_{i},w_{i}\rangle_{L^{2}(\nu)}=\int_{0}^{\tau}\omega_{i}(F(x))\omega_{j}(F(x))d\nu(x)$, and then we set $\tilde{\omega}_{i}(\cdot)=\sum_{j}\boldsymbol{P}^{-1}_{ij}w_{j}(\cdot)$. Notice that the matrix $\boldsymbol{P}$ may not have an inverse, meaning that the functions $\omega_{i}$ are not linearly independent in $L^{2}(\nu)$. In such a case we compute its Moore-Penrose pseudo-inverse. Since, in practice, we do not have enough information to compute $\boldsymbol{P}$, we can estimate it (under the null hypothesis) by $\widehat{\boldsymbol{P}}_{ij}=n/(n_{0}n_{1})\int_{0}^{\tau}\omega_{i}(\widehat{F}(x-))\omega_{j}(\widehat{F}(x-))L(x)d\widehat{\Lambda}(x)$. Finally, observe that this procedure generates a random kernel that converges to a deterministic one when $n$ tends to infinity.

4 Asymptotic Properties

In this section we study some asymptotic properties of our test statistic $Z_{n}$. Before introducing our main results, we state some standard and very reasonable conditions that are needed in our proofs. We first state some smoothness, and moment conditions regarding the kernel $K$.

Condition 3**.**

Let $K:[0,1)^{2}\to\mathds{R}$ be a continuous kernel, and let $X,X^{\prime},Y,Y^{\prime}$ be independent random variables such that $X,X^{\prime}\sim F_{0}$ and $Y,Y^{\prime}\sim F_{1}$. Then, we assume that

1. i)

$\mathds{E}\left(K(F(X),F(Y))^{2}(1-G_{0}(X))(1-G_{1}(Y))\right)<\infty$, 2. ii)

$\mathds{E}\left(K(F(X),F(X^{\prime}))^{2}(1-G_{0}(X))(1-G_{0}(X^{\prime}))\right)<\infty$, 3. iii)

$\mathds{E}\left(K(F(X),F(X))(1-G_{0}(X))\right)<\infty$, 4. iv)

$\mathds{E}\left(K(F(Y),F({Y^{\prime}}))^{2}(1-G_{1}(Y))(1-G_{1}(Y^{\prime}))\right)<\infty$, and 5. v)

$\mathds{E}\left(K(F(Y),F(Y))(1-G_{1}(Y))\right)<\infty$.

In the previous conditions, recall that $F=\eta F_{0}+(1-\eta)F_{1}$ is the pooled distribution, and $F=F_{0}$ under the null hypothesis. In addition, we need a technical condition to deal with the randomness generated by $\widehat{F}(x-)$, which is an analogue of the conditions needed in Theorem 4.2.1 of Gill (1980).

Condition 4**.**

Assume that for every $\varepsilon>0$, it holds that

1. i)

$\lim_{t\to\tau}\limsup_{n\to\infty}\mathds{P}\left(\int_{t}^{\tau}K(\widehat{F}(x-),\widehat{F}(x-))(1-G_{0}(x))dF_{0}(x)>\varepsilon\right)=0$, and 2. ii)

$\lim_{t\to\tau}\limsup_{n\to\infty}\mathds{P}\left(\int_{t}^{\tau}K(\widehat{F}(x-),\widehat{F}(x-))(1-G_{1}(x))dF_{1}(x)>\varepsilon\right)=0$.

Note that our conditions are not easy to verify in practice as they require knowledge of the distributions involved, however, the conditions are trivially satisfied for continuous and bounded kernels on $[0,1)^{2}$, which is a property that most standard kernels, such as the Gaussian kernel, satisfy.

Finally, in order to ease notation, we present our results in terms of the functions $\psi:\mathds{R}_{+}\to\mathds{R}$ and $\psi^{*}:\mathds{R}_{+}\to\mathds{R}$, defined by

[TABLE]

and

[TABLE]

4.1 Limit distribution under the null hypothesis

Our first result establishes that, under the null hypothesis, $(n_{0}n_{1}/n)Z_{n}$ converges in distribution to a limit random variable $Z$, when the number of data points tends to infinity.

Theorem 5** (Limit distribution under the null).**

Assume Conditions 3 and 4 hold. Then, under the null hypothesis, we have that

[TABLE]

as $n$ approaches infinity, where $Y=\sum_{i=1}^{\infty}\lambda_{i}(\xi_{i}^{2}-1)$, $\xi_{1},\xi_{2},\ldots$ are a collection of (potentially infinity) i.i.d. standard normal random variables, and $\lambda_{1},\lambda_{2},\ldots$ are non-negative constants. Additionally, the mean and variance of $Z$ are respectively given by

[TABLE]

In the previous theorem, the constants $\lambda_{i}$ are the eigenvalues of an integral operator in the $L^{2}$ space associated with the measure induced by the triple $(T_{i},\Delta_{i},c_{i})$ under the random covariates model. More details about the limit distribution are shown in Section B along with the proof of Theorem 5.

4.2 Power under alternatives

We continue by studying the asymptotic behaviour of our test statistic under the alternative hypothesis, that is, under the assumption that $H_{1}:F_{0}\neq F_{1}$ holds. To this end, we first prove that $Z_{n}$ (without considering any scaling) converges to a deterministic value.

Define the measures $\nu_{0}$ and $\nu_{1}$ on $\mathds{R}_{+}$ by

[TABLE]

where $\psi^{*}$ is defined in Equation (15), and define their embeddings into $\mathcal{H}$ by

[TABLE]

respectively. The measures $\nu_{0}$ and $\nu_{1}$ can be understood as the population measures, respectively, of the empirical measures $\nu_{0}^{n}$ and $\nu_{1}^{n}$ defined in Equation (11). We will show that, under appropriate conditions, the embeddings $\phi_{0}^{n}$ and $\phi_{1}^{n}$, defined in Equation (10), converge to $\phi_{0}$ and $\phi_{1}$, respectively, where convergence is with respect to the norm of $\mathcal{H}$. This fact together with Theorem 1 yields the following result:

Theorem 6**.**

Assume Conditions 3 and 4 hold, then

[TABLE]

At this point, it should be clear that our test is consistent under the alternative hypothesis if $\phi_{0}\neq\phi_{1}$ whenever $\nu_{0}\neq\nu_{1}$. We will show that if $F_{0}\neq F_{1}$ then $\nu_{0}\neq\nu_{1}$, however, this result does not immediately extrapolate to $\phi_{0}$ and $\phi_{1}$. A sufficient condition to ensure that $\phi_{0}$ and $\phi_{1}$ are different when $F_{0}\neq F_{1}$ is that the kernel $K$ is c-universal (see definition in Section 3.1).

Corollary 7** (Consistency).**

Suppose that Conditions 3 and 4 hold, and that the kernel $K$ is c-universal. Additionally, assume that $1-G_{c}(x)=0$ implies $S_{c}(x)=0$, for any $c\in\{0,1\}$. Then, under the alternative hypothesis, we have that $n_{0}n_{1}/nZ_{n}\to\infty$, which implies that

[TABLE]

where $Q_{n}(1-\alpha)$ denotes the $(1-\alpha)$-quantile of the random variable $(n_{0}n_{1}/n)Z_{n}$ under the null hypothesis.

We remark that the previous result uses the $(1-\alpha)$-quantile of $(n_{0}n_{1}/n)Z_{n}$ under the null hypothesis, which converges to a finite quantity due to Theorem 5.

Notice that Corollary 7 establishes consistency for all the alternatives under the assumption of a $c$-universal kernel $K$, nevertheless, even if the kernel $K$ is not $c$-universal, we can still ensure consistency for particular alternatives, as we show in the following proposition.

Proposition 8**.**

Assume that Conditions 3 and 4 hold, and suppose that $d\Lambda_{1}(x)=e^{\theta\bar{\omega}(F_{0}(x))}d\Lambda_{0}(x)$ for some $\bar{\omega}\in\mathcal{H}$ and $\theta\neq 0$. Then $(n_{0}n_{1}/n)Z_{n}\to\infty$, and thus the test is consistent for such an alternative.

Note that if $d\Lambda_{1}(x)=e^{\theta\bar{\omega}(F_{0}(x))}d\Lambda_{0}(x)$ for some $\bar{\omega}\in\mathcal{H}$ and $\theta\neq 0$, then $\liminf_{n\to\infty}\operatorname{LR}_{n}(\bar{\omega})^{2}>0$. In this case,

[TABLE]

Therefore $(n_{0}n_{1}/n)Z_{n}\to\infty$, and thus the test is consistent for such a particular alternative. In general, any argument that ensures consistency of $\operatorname{LR}_{n}(\omega)^{2}$ for some $\bar{\omega}$ in $\mathcal{H}$ also applies to our test statistic $Z_{n}$, as long as the limit distribution of $(n_{0}n_{1}/n)Z_{n}$ exists under the null hypothesis.

5 Wild Bootstrap Implementation

Given $\alpha\in(0,1)$, we construct a statistical test of level $\alpha$ by rejecting the null hypothesis if $(n_{0}n_{1}/n)Z_{n}>Q_{n}(1-\alpha)$, where $Q_{n}(1-\alpha)$ is the $(1-\alpha)$-quantile of the distribution of $(n_{0}n_{1}/n)Z_{n}$ under the null hypothesis. As the distribution of $(n_{0}n_{1}/n)Z_{n}$ is usually unknown, we can use its asymptotic distribution, given in Theorem 5, to approximate the rejection region. Unfortunately, excluding exceptional cases, the limiting distribution is rather complex, and thus computing the asymptotic $(1-\alpha)$- quantile is very hard. In order to carry out our test, we introduce a simple Wild Bootstrap implementation of our testing procedure.

Recall, from Equation (8), that our test statistic $Z_{n}$ can be written as follows:

[TABLE]

Notice that the expression given for $Z_{n}$ looks like a V-statistic. However, what prevents $Z_{n}$ from being a V-statistic is that it depends on the functions $\widehat{F}(x-)$, $Y_{0}(x)/n_{0}$ and $Y_{1}(x)/n_{1}$, and $n/(n_{0}n_{1})L(x)$, which are all random functions that depend on all data points. In Appendix B.1 and Appendix B.3, we prove that we can replace these functions by their respective limits plus some small additive error term that vanishes sufficiently fast as the number of observations grows. From this result, we deduce that $Z_{n}$ can be approximated by a V-statistic, suggesting that the standard Wild Bootstrap sampling scheme of Dehling and Mikosch (1994) can be used to approximate the asymptotic distribution of $n/(n_{0}n_{1})Z_{n}$.

Let $\mathcal{W}=(W_{i})_{i\in[n]}$ be a sequence of i.i.d. random variables satisfying that $\mathds{E}(W_{i})=0$ and $\mathds{E}(W_{i}^{2})=1$, and such that they are independent of any other source of randomness, in particular, from the observations. Then, we define the Wild Bootstrap statistic associated with $Z_{n}$ by

[TABLE]

Similarly to the expression given in Equation (9), we can write $Z_{n}^{\mathcal{W}}$ as

[TABLE]

where $\mathbf{W}\in\mathds{R}^{n}\times\mathds{R}^{n}$ is the diagonal matrix whose entries are given by $\mathbf{W_{jj}}=W_{j}$. Also, recall that $\mathbf{V}\in\mathds{R}^{n}$ is defined by $\mathbf{V}_{j}=L(X_{j})(-1)^{c_{j}}\Delta_{j}/Y_{c_{j}}(X_{j})$ and $\widehat{\mathbf{K}}\in\mathds{R}^{n}\times\mathds{R}^{n}$ is defined by $\widehat{\mathbf{K}}_{ij}=K(\widehat{F}(X_{i}),\widehat{F}(X_{j}))$. This matrix expression is very convenient for the computational implementation of our method.

The following theorem states the correctness of our Wild Bootstrap approach:

Theorem 9**.**

Assume that Conditions 3 and 4 are satisfied. Then, under the null hypothesis, it holds that for any $x\in\mathds{R}_{+}$,

[TABLE]

where $Z$ is the random variable defined in Theorem 5.

With the previous ingredients we are ready to describe our testing procedure, which relies on approximating the $(1-\alpha)$-quantile of the distribution of $Z_{n}$ by the quantiles of $Z_{n}^{\mathcal{W}}$. Since we can freely sample independent copies from $Z_{n}^{\mathcal{W}}$ given the data points, we can estimate the quantile by Monte Carlo simulations. Our algorithm is as follows:

1. i)

Set level $\alpha\in(0,1)$ of the test and let $N$ be a large integer, 2. ii)

Sample $N$ independent copies of the Wild Bootstrap statistic using Equation (22), 3. iii)

Compute the $(1-\alpha)$-quantile of the previous sample and call it $Q_{n}^{\mathcal{W}}(1-\alpha)$, 4. iv)

Compute the test statistic $Z_{n}$ using Equation (9), 5. v)

Reject the null-hypothesis if $Z_{n}>Q_{n}^{\mathcal{W}}(1-\alpha)$, otherwise do not reject it.

6 Simulations

We perform an empirical evaluation of our methods in which the ground truth is known. To this end, we consider two different settings: a proportional hazard functions setting (in which the classic log-rank test is provably the most powerful), and a time-dependent hazard functions setting, including Weibull and periodic hazard functions. All our experiments consider the same null cumulative hazard function $\Lambda_{0}(t)=t$, corresponding to the cumulative hazard function of an exponential random variable with mean 1. We choose $\Lambda_{1}(t)$ belonging to one of the following parametric families:

1. i)

Proportional hazards: for this case we consider $\Lambda_{1}(t)$ belonging to the parametric family given by $\Lambda(t;\theta)=\theta t$ with $\theta\in[0,2]$. Observe that $\Lambda_{0}(t)$ is recovered when $\theta=1$. 2. ii)

Weibull (polynomial) hazards: we consider $\Lambda_{1}(t)$ belonging to a family of Weibull cumulative hazard functions $\Lambda(t;\theta)=t^{\theta}$ with $\theta\in[0,2]$. Notice that $\theta=1$ recovers the null. 3. iii)

Periodic hazards: we consider $\Lambda_{1}(t)$ belonging to the family given by $\Lambda(t;\theta)=t-\sin(\pi\theta t)/(\pi\theta)$ with $\theta\in(0,15]$. Notice that $\lim_{\theta\to\infty}\Lambda(t;\theta)=\Lambda_{0}(t)=t$, recovering the null hypothesis.

Figure 1 shows the behaviour of the cumulative hazard functions of the parametric families previously described, for different values of $\theta$.

6.1 Implementation

6.1.1 Kernels

The heart of our testing-approach is undoubtedly the kernel function. In Section 3.3 we showed that specific choices of the kernel function lead to some existing tests. In our experiments we use some of these choices as well as a kernel that is $c$-universal. We define kernels in the following categories:

Log-rank kernels (LRP and LRC): In Section 3.3.1, we considered kernels of the form $K(x,y)=\omega(x)\omega(y)$, recovering the well-known weighted log-rank tests. In our experiment we choose $\omega$ to be equal to i) $\omega_{1}(x)=1$ and ii) $\omega_{2}(x)=(x-1/2)$. For i) we recover the classical log-rank test (LRP), used to test proportional hazard functions, while for ii) we recover a weighted log-rank test (LRC), designed to discover if the hazard functions cross around the median of the distribution $F$. 2. 2.

Projection kernels (P2W and P4W): we follow the approach described in Section 3.3.3, which recovers the testing procedure of Brendel et al. (2014). In particular, we choose kernel functions based on the subspace generated by the weight functions i) $\{1,x\}$ and ii) $\{1,x,x^{2},x^{3}\}$. We denote the kernels in i) and ii) by P2W and P4W, respectively, making a clear reference to the dimension of the subspace. 3. 3.

Pearson-type kernels (Per4 and Per5): we consider the Pearson-type kernel functions defined in Section 3.3.2, particularly in Equation (13). For our experiments, we partition the space into 4 and 5 disjoint regions, and denote the kernels by Per4 and Per5, respectively. 4. 4.

Squared exponential kernel (SEK): we consider the squared exponential kernel (SEK) defined by $K(x,y)=\exp\{-(x-y)^{2}/\sigma^{2}\}$ with $\sigma=0.1$. This kernel is $c$-universal, hence by Corollary 7, it leads to an omnibus test. Better results may be obtained by optimising the parameter $\sigma$ for the problem at hand. A well-known heuristic is to choose $\sigma$ as the median of the pairwise differences of the observations Scholkopf and Smola (2001).

.

6.1.2 Computer Implementation

Our experiments are implemented in R following the Wild Bootstrap approach described in Section 5, choosing Rademacher random variables for the Wild Bootstrap weights $\mathcal{W}$. We compute our test statistic $Z_{n}$ and the Wild Bootstrap statistic $Z_{n}^{\mathcal{W}}$ by using the quadratic form expressions given in Equations (9) and (22), respectively. By using these expressions, it is possible to give a simple and fast implementation of our testing procedure, indeed, 1000 repetitions of our testing procedure (using 1000 Wild Bootstrap samples) takes just a couple of minutes in a standard commercial laptop, for $n_{0}=n_{1}=100$. For each experiment, we consider sample sizes of $30$ and $100$ observations per group, and we choose a censoring distribution generating $10\%$ and $30\%$ of censored observations.

6.2 Type-I error

In our first experiment, we verify that our tests achieve a correct Type-I error of $\alpha=0.05$ for each combination of sample size and censoring percentage, for the kernels described in Section 6.1.1. For each different combination of parameters, we run our test in 1000 simulated datasets. Table 1 shows the results. In general, the Wild Bootstrap approach has no problem reaching the correct level for even censoring percentages. Uneven censoring causes a few problems, so the user should be careful when applying this method in the latter setting. Arguably, the squared exponential kernel (SEK) is the most robust. It is worth recalling that projection tests use a random kernel (it depends on the data points) which might impact its performance. Additionally, while not reported here, we tested other significance levels, obtaining similar results.

6.3 Power Simulations

We provide an empirical evaluation of our testing procedure for each of the settings described in Section 6: proportional, Weibull and periodic hazard functions. As previously mentioned, the null cumulative hazard function is given by $\Lambda_{0}(t)=t$. The power is estimated by repeating our testing procedure over 1000 simulated datasets for each combination of sample size and censoring percentage. Results for the proportional, Weibull and periodic settings are shown in Figures 2, 3, and 4, respectively. We give a few comments and remarks about our experiments.

We report a very small fraction of our experiments as all the results are qualitatively the same as the ones shown in Figures 2, 3, and 4. 2. 2.

We only report results for the kernel Per5 since the kernel Per4 has an almost identical behaviour. 3. 3.

The LRP test is equivalent to the score test for the proportional hazard functions model; thus, it can be deduced that it is the most powerful test for local alternatives under this model. The previous statement is supported by Figure 2, where we observe an excellent performance of the LRP test in the setting of proportional hazard functions alternatives. Observe that the LRP test loses all of its power (nearly zero power) for the Weibull and periodic hazard alternatives as shown in Figures 3, and 4. 4. 4.

The LRC test is equivalent to the weighted log-rank test with weight $w(x)=x-1/2$. It can be observed in Figure 3, that the LRC test has a very good performance in the setting of Weibull hazard alternatives, but its power is relatively low in the other two settings, as shown in Figures 2 and 4. An explanation for this behaviour follows from the fact that the LRC kernel is designed to be optimal at detecting a cross around the median of the pooled distribution $F$, and thus, it will not be a good fit for the proportional hazards experiment. Also, for more complex hazard functions, such as those described in the periodic setting, we can observe more than one cross occurring, which explains the poor performance of the LRC test. 5. 5.

We observe different behaviours for the projection kernels P2W and P4W. In the proportional hazard functions setting, the P2W test has the best performance after the LRP test, which is explained by the fact that this kernel is constructed considering the subspace generated by $\{1,x\}$, where the weight $\omega(x)=1$ is known to be optimal for proportional hazard alternatives. While P4W also includes this weight (it is generated by $\{1,x,x^{2},x^{3}\}$), the fact that it considers a larger space of possible alternatives makes the test more data-expensive resulting in a loss of power. For the Weibull hazard functions setting, we observe that both projection tests, P2W and P4W, have an overall good performance. This behaviour can be explained by the fact that both tests consider projections on polynomials, and the Weibull hazard functions are, indeed, polynomials. In the periodic case, both kernels have a substandard behaviour due to the fact that the hazard structure is rather different to a polynomial of finite degree. Note that with more data it seems that the tests do not improve when compared to the best kernel, in this case, the squared exponential kernel SEK. 6. 6.

The Pearson-type kernel Per5 has a consistent behaviour, being neither too good nor too bad. Disadvantages are that the user needs to specify beforehand a partition of the space. 7. 7.

The squared exponential kernel, SEK, gives overall good results. In one hand, while in the presence of more ‘structured’ data, i.e., proportional or Weibull (polynomial) hazard functions, simpler kernels give better results, the SEK still has a good performance. On the other hand, in presence of more ‘complex’ data, its performance is better than other kernels. In general, the SEK performs better than Pearson-type kernels, suggesting that the SEK is a better choice for an all-around kernel. 8. 8.

In general, it seems that for nice structured data, simpler kernels (leading to simpler methods) have better results. On the other hand, for complex data, a more complex kernel seems to be a good option. Overall, we think that the SEK is the best option as it has a robust behaviour in simple data-scenarios, and it outperforms other tests in more complex scenarios. Also, as shown in Table 1, this kernel has close to no problems reaching the correct Type-I error.

7 Real data

We consider the Gastrointestinal Tumor Study Group data (GTSG), Stablein et al. (1981), available in the ‘coin’ R-package. The data considers a randomised clinical trial in the treatment of locally advanced, non-resectable gastric carcinoma. In this study, 42 patients were treated by using chemotherapy alone, while 45 patients were treated by using a combination of chemotherapy and radiation therapy. The aim of the study is to detect differences between these treatments. Kaplan-Meier curves for each group are shown in Figure 5. The null hypothesis is that there is no difference between the treatments. We apply our test considering the 7 different kernels described in Section 6.1.1. The corresponding p-values (approximated by using our Wild Bootstrap approach) are shown in Table 2.

All the tests, except for the classical log-rank test, reject the null hypothesis at $5\%$ level of significance. This outcome is quite reasonable as the survival functions cross, as shown in Figure 5. Indeed, the smallest $p$-values are given by the SEK, LRC and the P2W tests which reject the null hypothesis at $1\%$ level of significance. This is not a surprising behaviour of LRC and P2W, as these kernels are tailored to detect crossings. The SEK also performs very well which is quite satisfying as this kernel is not particularly designed for the setting of crossing hazards.

8 Conclusion

We have introduced an RKHS-based testing procedure for the standard problem of two-sample hypothesis testing in the framework of right censored data. Our test statistic is the supremum of weighted log-rank statistics, where the weights belong to the unit ball of an RHKS. While our test statistic is apparently very complex, its evaluation becomes analytically tractable due to the reproducing property of reproducing kernel Hilbert spaces. Indeed, our test statistic can be written as a quadratic form as in Equation (9), and it can be fully characterised by a kernel function $K$. This simple structure allows us to derive asymptotic properties of our test, and it suggests that the standard Wild bootstrap approach can be used to approximate the rejection region. We also showed that our test can be seen as a natural infinitely-dimensional generalisation of other well-known test statistics based on log-rank statistics. Finally, we performed a simulation study which compares the results of the test statistic for different kernel functions.

To finish the paper we discuss some of our results and potential research ideas. First, as shown in Theorem (5), the asymptotic distribution of our test statistic under the null hypothesis is, in general, intractable, and thus we need to rely on bootstrapping techniques to approximate the rejection region. While, for us, the natural option is to consider a Wild Bootstrap approach, it is also possible to consider a permutation approach as the one used for weighted log-rank tests (Neuhaus, 1993). In the case of equal censoring, the permutation approach has the advantage of being finitely exact, however, in the general case of censoring, the permutation approach is not directly applicable since the limit distribution of our test statistic depends on the censoring distribution of each group. For weighted log-rank tests, this problem is fixed by standardising the test statistic leading to an asymptotically distribution-free test. Unfortunately, it is not clear how to standardise (if possible) our test statistic, and we think this is a non-trivial task, especially for infinitely-dimensional kernels, hence, we leave this problem as future research.

A natural question to ask is: How do we choose a kernel function? Unfortunately, nobody can provide an answer for such a question yet and choosing kernel functions is an active field of research in Statistics and Machine Learning, as it is an issue that happens in several other contexts such as Gaussian Processes inference, support vector machine, kernel regression, etc. In practice, our simulation study suggests that simple (finite-dimensional) kernels, e.g., polynomials, are good enough if we know in advance that the hazard function has a relatively simple form, whereas complex kernels are better suited for complex hazard functions. In general, the squared exponential kernel is a very safe choice as it performs well in both settings. Nevertheless, we should not expect it to perform well in every setting, as Janssen (2000) proved that for a finite number of data points, any nonparametric test has preferences for a finite-dimensional subspace of alternatives, and outside this subspace, the power of the test is almost flat in balls of alternatives. This means that, for a fixed number of data points, it is possible to find alternatives for which the squared exponential kernel has fails, however, due to the complex form of the kernel, it may be difficult to construct such an alternative, and we do not expect that to happen in practical settings. Note that, by choosing a reproducing kernel, Theorem 5 shows that our test fixes most of its power on a finite number of directions, which are given by the eigenfunctions with larger eigenvalues in the spectral decomposition of the kernel, and thus, with a finite number of data points we expect the test to concentrate its power around the first $m$ eigenfunctions with larger eigenvalues, where $m$ is some constant depending on $n$, and as long as $n$ grows, $m$ should be growing as well.

Finally, we give a few comments about our technical results. Our asymptotic analysis uses the fact that our test statistic is a double stochastic integral (Equation (7)), and particularly, under the null hypothesis, those stochastic integrals are with respect to martingales. Our main tool to study these objects is Theorem 17, which allows us to control double stochastic integrals with respect to a special class of random integrands $h$, which includes our integrand $K(\widehat{F}(x-),\widehat{F}(y-))L(x)L(y)$, as well as several others. Our theorem is powerful enough to study other type of collection of log-rank statistics, for example, we can consider the test statistic

[TABLE]

which is similar to our test statistic (c.f., Equations (4) and (1)), but replacing $\widehat{S}(x-)$ by $Y(x)/n$ (See (Fleming and Harrington, 1991, Chapter 7)). In general, replacing $\widehat{S}(x-)$ by any reasonable predictable process gives a test statistic that can be analysed by our methods, and furthermore, our tool can be applicable to an even larger class of processes. The analysis of other type of test statistics is subject of future research.

Correspondence

Nicolas Rivera.

Address: Computer Laboratory, University of Cambridge. William Gates Building, 15 JJ Thomson Ave, Cambridge CB3 0FD, UK.

e-mail: nr454 at cam dot ac dot uk

Appendix

Appendix A Preliminary Results

A.1 Some Results for Counting Processes

Proposition 10**.**

The following results hold:

1. i)

$\lim_{n\to\infty}\sup_{t\leq\tau}|\widehat{S}(t)-S(t)|=0$* a.s.,* 2. ii)

$\lim_{n\to\infty}\sup_{t\leq\tau}|Y(t)/n-H(t-)|=0$* a.s.,*

Item i is due to Stute and Wang [1993], and item ii is the Glivenko-Cantelli theorem.

Proposition 11**.**

Let $\beta\in(0,1)$, then

1. i)

$\mathds{P}(\widehat{S}(t)\leq\beta^{-1}S(t),\quad\forall t\leq\tau_{n})\geq 1-\beta,$** 2. ii)

$\mathds{P}\left(Y(t)/n\leq\beta^{-1}(1-H(t-)),\quad\forall t\leq\tau_{n}\right)\geq 1-\beta,$* and* 3. iii)

$\mathds{P}\left(Y(t)/n\geq\beta(1-H(t-)),\quad\forall t\leq\tau_{n}\right)\geq 1-e(1/\beta)e^{-1/\beta}.$**

i.e., i) $\sup_{t\leq\tau_{n}}\widehat{S}(t)/S(t)=O_{p}(1)$, ii) and iii) $\sup_{t\leq\tau_{n}}Y(t)/(n(1-H(t-))=\Theta_{p}(1)$.

The proofs of items i and iii are due Gill et al. [1983], and item ii follows from Gill [1980, Theorem 3.2.1].

Remark 12*.*

Notice that the results of the previous propositions still hold if we replace $Y(t)$, $H(t)$ and $n$ by $Y_{c}(t),H_{c}(t)$ and $n_{c}$, respectively.

A.2 Double Martingales

In our proofs, we will frequently encounter double martingale integral processes of the form:

[TABLE]

where $h_{n}:\mathds{R}^{2}\to\mathds{R}$ is a sequence of symmetric positive-definite random functions, and $(W_{n}(t))_{t\geq 0}$ is a sequence of $(\mathcal{F}_{t})$-martingales. The aim of this section is to establish conditions under which such processes converge to zero in probability when evaluated at $t=\tau_{n}$. This result is formally stated in Theorem 17, and to prove it, we use some results introduced by Fernández and Rivera [2020] for double integrals with respect to martingales.

Definition 13**.**

Define the $\sigma$-algebra $\mathcal{P}$ on $C=\{(x,y):0<x<y<\infty\}$ as the $\sigma$-algebra generated by sets of the form

[TABLE]

and $\{0\}\times\{0\}\times X$, where $X\in\mathcal{F}_{0}$.

A stochastic process $(h(x,y))_{(x,y)\in C}$ is said to be $\mathcal{P}$-measurable if it is measurable with respect to $\mathcal{P}$.

The following proposition is a simple consequence of the definition of $\mathcal{P}$.

Proposition 14**.**

Let $k:\mathds{R}^{2}\to\mathds{R}$ be a measurable function, and let $(h_{1}(t))_{t\geq 0}$ and $(h_{2}(t))_{t\geq 0}$ be $(\mathcal{F}_{t})$-predictable stochastic processes. Then the process $(h(x,y))_{(x,y)\in C}$ given by $h(x,y)=k(h_{1}(x),h_{2}(y)))$ is $\mathcal{P}$-measurable.

Theorem 15**.**

Let $h$ be a $\mathcal{P}$-measurable process, and let $W$ be a right-continuous $(\mathcal{F}_{t})$-martingale. Assume that for all $t\geq 0$, it holds that

[TABLE]

Then, $(Z(t))_{t\geq 0}$, defined by $Z(t)=\int_{0}^{t}\int_{(0,y)}h(x,y)dW(x)dW(y)$, is an $(\mathcal{F}_{t})$-martingale.

In our proofs we are particularly interested in predictable positive-definite processes which are defined as following:

Definition 16**.**

We say a process $(h(x,y))_{x,y\geq 0}$ is a predictable positive-definite process if it satisfies the following properties: i) $h(x,y)=h(y,x)$, ii) $h$ is positive definite, i.e., each realisation of the stochastic process is a positive definite function, iii) $(h(x,y))_{(x,y)\in C}$ is $\mathcal{P}$-measurable, and iv) $(h(x,x))_{x\geq 0}$ is predictable with respect to $(\mathcal{F}_{x})_{x\geq 0}$.

The next theorem, whose proof is deferred to Section E, gives sufficient conditions under which the process of Equation (23) converges to zero in probability.

Theorem 17**.**

Let $(h_{n}(x,y))_{x,y\geq 0}$ be a sequence of predictable positive-definite processes, let $W_{n}$ be a sequence of right-continuous $(\mathcal{F}_{t})$-martingales with predictable and quadratic variation processes denoted by $\langle W_{n}\rangle$ and $[W_{n}]$, respectively, and suppose that

[TABLE]

holds for all $n$ large enough, and for all $t\in(0,\tau)$.

Then, if

[TABLE]

we have that

[TABLE]

and the same holds if we replace $o_{p}(1)$ by $O_{p}(1)$.

Note that Equation (25) holds trivially due to the simple nature of our martingales, hence, we will not verify this conditions in our applications of the theorem.

Appendix B Analysis under the Null Hypothesis: Proof of Theorem 5

The proof of Theorem 5 is split into three mains steps:

• i)

We find a cleaner asymptotic expression for our test statistic under the null hypothesis. In particular, we show that $Z_{n}$ can be rewritten as

[TABLE]

Then, by using Theorem 17, we prove that $\widehat{F}(x-)$, $n/(n_{0}n_{1})L(x)$, and $Y_{c}(x)/n$ can be replaced by their respective limits, $F_{0}(x)$, $\psi(x)S_{0}(x)$, and $(1-H_{c})(x)\eta^{c-1}(1-\eta)^{-c}$, up to a small additive term that decreases to zero in probability, obtaining that

[TABLE]

• ii)

We prove that the deterministic and random covariates models are asymptotically equivalent in the sense that our test statistic $nZ_{n}$ has the same asymptotic distribution (when it exists) under both models.

• iii)

We obtain the limit distribution of $nZ_{n}$ under the random covariates model. Our results translate to the deterministic covariates model by using the result of the previous item.

B.1 Step i: Finding a simpler asymptotic representation

For this step, we work under the deterministic covariates model, but notice that our analysis can be extended to the random covariates model by conditioning on the number of random covariates with value 0, say $N_{0}$, and by noticing that $N_{0}/n\to\eta$ almost surely.

From Theorem 1 and Lemma 2, we have that

[TABLE]

where $d\Lambda^{*}(x)=d\widehat{\Lambda}_{0}(x)-d\widehat{\Lambda}_{1}(x)$, and

[TABLE]

holds under the null hypothesis. In particular, notice that, under the null hypothesis, $(\Lambda^{*}(x))_{x\geq 0}$ is an $(\mathcal{F}_{x})$-martingale with predictable and quadratic variation processes given by

[TABLE]

respectively. Then, by substituting Equation (29) in Equation (28), we can rewrite our test statistic as

[TABLE]

The main result of this section is the following:

Theorem 18**.**

Assume Conditions 3 and 4 hold. Then, under the null hypothesis, it holds that

[TABLE]

We split the proof of Theorem 18 into three parts:

We prove that the pooled Kaplan-Meier estimator, $\widehat{F}(x-)$ can be replaced by its limit $F_{0}(x)$ (recall $F_{0}$ is continuous, and thus $F_{0}(x)=F_{0}(x-))$, 2. 2)

we prove that $nL(x)/(n_{0}n_{1})$ can be replaced by its limit $\psi(x)S_{0}(x)$, and 3. 3)

we prove that $Y_{0}(x)/n$ and $Y_{1}(x)/n$ can be replaced by their limits, $\eta(1-H_{0}(x))$ and $(1-\eta)(1-H_{1}(x))$, respectively.

For the proof of Theorem 18 we will recurrently use the following fact: let $q:\mathds{R}\to[0,1)$ and let $\mu$ be a measure on $[0,\infty)$, then

[TABLE]

which is a straightforward consequence of the linearity of the inner product and integration, and the reproducing property (this argument was used in the proof of Theorem 1).

B.2 Proof of Theorem 18

**Part 1): Replacement of $\widehat{F}(x-)$ by $F_{0}(x)$

**The following Lemma proves that $K(\widehat{F}(x-),\widehat{F}(y-))$ can be replaced by $K(F_{0}(x),F_{0}(y))$, up to a small error.

Lemma 19**.**

Assume Conditions 3 and 4. Then, under the null hypothesis, it holds that

[TABLE]

Proof.

Using norm notation, the desired result is equivalent to show that

[TABLE]

By triangular inequality $\|b\|_{\mathcal{H}}-\|a-b\|_{\mathcal{H}}\leq\|a\|_{\mathcal{H}}\leq\|b\|_{\mathcal{H}}+\|a-b\|_{\mathcal{H}}$, then, by taking square, we deduce that we just need to prove that

[TABLE]

and that

[TABLE]

We begin by verifying Equation (33). Expanding the inner product expression, the left-hand side equals

[TABLE]

Also, notice that $(K(F(x),F(y))L(x)L(y))_{x,y\geq 0}$ is a predictable positive-definite process (recall Definition 16), and since $(\Lambda^{*}(x))_{x\geq 0}$ is an $(\mathcal{F}_{x})$-martingale with predictable and quadratic variation processes given in Equation (30), a straightforward application of Theorem 17 tell us that we just need to verify that

[TABLE]

The previous equation holds true by Condition 3, and by using that $L(x)/n=O_{p}(1)S_{0}(x)\psi(x)$ uniformly for all $x\leq\tau_{n}$ due to Propositions 11.ii and 11.iii, and since $\psi(x)\leq\eta^{-1}(1-G_{0}(x))$.

We continue verifying Equation (32). Notice that

[TABLE]

where

[TABLE]

It is easy to verify that $(h(x,y))_{x,y\geq 0}$ is a predictable positive-definite process, then, by using Theorem 17, the desired result follows from proving $n\left(\frac{n}{n_{0}n_{1}}\right)^{2}\int_{0}^{\tau_{n}}h(x,x)d\langle\Lambda^{*}\rangle(x)=o_{p}(1)$.

Using that $L(x)/n=O_{p}(1)S_{0}(x)\psi(x)$ uniformly for all $x\leq\tau_{n}$ and $\psi(x)\leq\eta^{-1}(1-G_{0}(x))$, we get

[TABLE]

where $\omega(x)=K(F_{0}(x),F_{0}(x))-2K(F_{0}(x),\widehat{F}(x-))+K(\widehat{F}(x-),\widehat{F}(x-)))$. Let $\varepsilon>0$, then for any $t\in(0,\tau)$,

[TABLE]

We will prove that both terms on the right-hand side of the previous equation tends to 0 as $n$ approaches infinity. For the first term, notice that, since $\sup_{x\leq t}|\widehat{F}(x-)-F(x)|=0$ a.s. by Proposition 10.i, and since $K$ is continuous in $[0,1)^{2}$, it exists $N$, large enough such that for every $n\geq N$, $\omega(x)\leq\varepsilon/4$ uniformly on $[0,t]$, thus $\int_{0}^{t}\omega(x)(1-G_{0}(x))dF_{0}(x)\leq\varepsilon/4$. Therefore, for any $t<\tau$,

[TABLE]

For the second term of the right-hand side of Equation (35), note that $2K(F_{0}(x),\widehat{F}(x-))\leq K(F_{0}(x),F_{0}(x))+K(\widehat{F}(x-)),\widehat{F}(x-))$, then $\omega(x)\leq 2K(F_{0}(x),F_{0}(x))+2K(\widehat{F}(x-)),\widehat{F}(x-))$, and thus we deduce

[TABLE]

By Condition 3 we deduce that it exists $t^{\prime}<\tau$, $\int_{t^{\prime}}^{\tau}K(F_{0}(x),F_{0}(x)(1-G_{0}(x))dF_{0}(x)<\varepsilon/16$. Then, by combining Equations (36) and (37), we get that for every $t\geq t^{\prime}$

[TABLE]

and by taking $t\to\tau$, we deduce that $\limsup_{n\to\infty}\mathds{P}\left(\int_{0}^{\tau}\omega(x)(1-G_{0}(x))dF_{0}(x)\geq\varepsilon\right)=0$ by Condition (4), completing our proof. ∎

**Part 2): Replacement of $\frac{n}{n_{0}n_{1}}L(x)$ by $\psi(x)S_{0}(x)$

**

By the previous part, our test statistic satisfies

[TABLE]

The next step is to show that we can replace $n/(n_{0}n_{1})L(x)$ by its limit $\psi(x)S_{0}(x)$ without altering the value of $nZ_{n}$ by more than an error of order $o_{p}(1)$. We formalise this result in the following lemma.

Lemma 20**.**

Under Conditions 3 and 4, it holds

[TABLE]

Proof.

Similar to the proof of Lemma 19, we just need to prove that

[TABLE]

and

[TABLE]

We begin by proving Equation (38). Note that, by expanding the inner product, it is enough to prove that

[TABLE]

where

[TABLE]

Also, observe that $(h(x,y))_{x,y\geq 0}$ is a predictable positive-definite process. Then, by Theorem 17, we just need to verify that $n\int_{0}^{\tau_{n}}h(x,y)d\langle\Lambda^{*}\rangle(x)=o_{p}(1),$ which is deduced from the following equalities:

[TABLE]

where the second equality is due to Propositions 11.ii and 11.iii, from which we deduce that $L(x)/n=\Theta_{p}(1)\psi(x)S_{0}(x)$ uniformly for all $x\leq\tau_{n}$, and the third equality is due to an application of the dominated convergence theorem in sets of arbitrarily large probability: observe that, by Proposition 10, we have that $\left(1-\frac{\psi(x)S_{0}(x)}{n/(n_{0}n_{1})L(x)}\right)\to 0$ for all $x<\tau$. Additionally,

[TABLE]

since $\psi(x)\leq(1-G_{0}(x))$ for all $x<\tau$, and note that $K(F(x),F(x))(1-G_{0}(x))$ is integrable by Condition 3. With these ingredients we can apply the dominated convergence theorem in set of arbitrarity large probability.

To check Equation (39), we follow the same steps, using that $n\int_{0}^{\tau_{n}}K(F(x),F(x))\psi(x)^{2}S_{0}(x)^{2}/L_{0}(x)d\Lambda_{0}(x)=O_{p}(1)$, since $L_{0}(x)/n=\Theta_{p}(1)\psi(x)S_{0}(x)$ uniformly for all $x\leq\tau_{n}$. ∎

**Part 3): Replacement of $Y_{0}$ by $\eta n(1-H_{0})$ and $Y_{1}$ by $(1-\eta)n(1-H_{1})$ ** From the previous step, it holds that

[TABLE]

Our next step is to replace $Y_{0}/n$ and $Y_{1}/n$ by their corresponding limits $\eta S_{0}(1-G_{0})$ and $(1-\eta)S_{0}(1-G_{1})$. The following Lemma implies the desired result.

Lemma 21**.**

Under Conditions 3, it holds

[TABLE]

and

[TABLE]

for any class label $c\in\{0,1\}$.

Proof.

We only prove the result for $c=0$ as the proof for $c=1$ follows the same steps. Define $\sigma_{n}(x)=\left(\frac{1}{Y_{0}(x)/n}-\frac{\eta^{-1}}{S_{0}(x)(1-G_{0}(x))}\right),$ and notice that

[TABLE]

Define $h(x,y)$ as the integrand in the previous double integral, and notice that $(h(x,y))_{x,y\geq 0}$ is a predictable positive-definite process. Also, recall that $M_{0}(x)$ is an $(\mathcal{F}_{x})_{x\geq 0}$-martingale with predictable and quadratic variation processes given by $d\langle M_{0}\rangle(x)=Y_{0}(x)d\Lambda_{0}(x)$ and $d[M_{0}](x)=dN(x)$. Then, by Theorem 17, we just need to verify that

[TABLE]

By Proposition 11.ii, we have $Y_{0}(x)/n=O_{p}(1)S_{0}(x)(1-G_{0}(x))$ uniformly for all $x\leq\tau_{n}$.

[TABLE]

To prove that the previous expression is $o_{p}(1)$ we use the dominated convergence theorem. For such, observe that $\sigma_{n}(x)\to 0$ a.s. for all $x<\tau$, and that $\sigma_{n}(x)^{2}=O_{p}(1)(S_{0}(x)(1-G_{0}(x))^{2}$ uniformly for all $x\leq\tau_{n}$ due to proposition 11.iii, therefore

[TABLE]

uniformly for all $x\leq\tau$. Finally, note that $K(F_{0}(x),F_{0}(x))\frac{\psi(x)^{2}}{1-G_{0}(x)}$ is integrable by Condition 3, since $\psi(x)\leq\min\{\eta^{-1}(1-G_{1}(x)),(1-\eta)^{-1}(1-G_{0}(x))\}$. The previous analysis yields that we can use the dominated convergence theorem in sets of large probability, concluding that Equation (40) holds true. Equation (41) follows from the same arguments, by noticing that $\int_{0}^{\tau_{n}}K(F_{0}(x),F_{0}(x))\frac{\psi(x)}{(1-G_{0}(x))}Y_{0}(x)d\Lambda_{0}(x)=O_{p}(1)$. ∎

B.3 Step II: Asymptotically equivalence of two models

After applying Theorem 18 we get

[TABLE]

Recall that the previous expression is valid when considering either deterministic or independent Bernoulli$(\eta)$ covariates. In this section we show that asymptotic results obtained using either the deterministic or random covariates models are equivalent.

Let $s_{i}=(X_{i},\Delta_{i},c_{i})$ for all $i\in[n]$ and define $J:(\mathds{R}_{+}\times\{0,1\}\times\{0,1\})^{2}\to\mathds{R}$ as

[TABLE]

Also, define

[TABLE]

and notice that $V_{n}$ is the non-negligible part of $Z_{n}$, that is, $Z_{n}=V_{n}+o_{p}(n^{-1})$ for both, the deterministic and the random covariates models. In this section only, we denote by $Z_{n}^{\prime}$ and $Z_{n}$, respectively, the test statistics under the deterministic and random covariates models. More generally, we use an apostrophe (e.g., $V_{n}^{\prime}$) to denote any term related to the deterministic covariates model.

Lemma 22**.**

There exists a coupling of $Z_{n}$ and $Z_{n}^{\prime}$ such that $Z_{n}-Z_{n}^{\prime}=o_{p}(n^{-1})$.

Proof.

We construct a coupling such that $V_{n}-V_{n}^{\prime}=o(n^{-1})$, which implies the desired result. The coupling is constructed as follows:

For $c\in\{0,1\}$, denote by $\mu_{c}$ the measure induced by the random pair $(X,\Delta)$ where $X=\min(T,C)$ and $\Delta=\mathds{1}_{\{T\leq C\}}$, and $T$ and $C$ are independent random variables such that $T\sim F_{c}$ and $C\sim G_{c}$. 2. 2.

Let $\mathcal{L}_{0}=(X_{i},\Delta_{i},0)_{i\geq 1}$ and $\mathcal{L}_{1}=(\bar{X}_{i},\bar{\Delta}_{i},1)_{i\geq 1}$ be such that all the triples in $\mathcal{L}_{0}$ and $\mathcal{L}_{1}$ are independent of everything, and $(X_{i},\Delta_{i})\sim\mu_{0}$ and $(\bar{X}_{i},\bar{\Delta}_{i})\sim\mu_{1}$ for all $i\geq 1$. Note that, since the times are sampled from continuous distributions, the data points are unique almost surely. 3. 3.

In the deterministic covariates model we generate the data, namely $\mathcal{D}^{\prime}_{n}$, by choosing the first $n_{0}$ elements of $\mathcal{L}_{0}$ and the first $n_{1}$ elements of $\mathcal{L}_{1}$. 4. 4.

In the random covariates model, we generate the data, namely $\mathcal{D}_{n}$, by sampling $N_{0}\sim\text{Binomial}(n,\eta)$ and then choosing the first $N_{0}$ elements of $\mathcal{L}_{0}$ and the first $n-N_{0}$ elements of $\mathcal{L}_{1}$. 5. 5.

Compute $V_{n}^{\prime}$ and $V_{n}$ by using the datasets $\mathcal{D}^{\prime}_{n}$ and $\mathcal{D}_{n}$, respectively.

Note that the datasets $\mathcal{D}_{n}$ and $\mathcal{D}_{n}^{\prime}$ differ in exactly $Q=|N_{0}-n_{0}|$ points. For ease of notation, denote by $s_{i}=(X_{i},\Delta_{i},c_{i})$ the observations in $\mathcal{D}_{n}$ and by $s_{i}^{\prime}=(X_{i}^{\prime},\Delta_{i}^{\prime},c_{i}^{\prime})$ the observations in $\mathcal{D}_{n}^{\prime}$. We sort the datasets in such a way that the first $Q$ elements of each dataset are not present in the other, and the remaining $n-Q$ elements in such a way that $s_{j}=s_{j}^{\prime}$ for $j\geq Q+1$. Then, a simple computation shows that

[TABLE]

where

[TABLE]

We continue by proving that $nD_{n}=o_{p}(1)$ and $nR_{n}=o_{p}(1)$. To this end, we need the following intermediate result, whose proof is deferred to Appendix E.

Proposition 23**.**

Let $A:\mathds{R}_{+}\times\mathds{R}_{+}\to\mathds{R}$ be a symmetric function such that, for every $c,c^{\prime}\in\{0,1\}$, the following integrals are finite: $\int_{0}^{\tau}|A(x,x)|(1-G_{c}(x))dF_{c}(x)$ and $\int_{0}^{\tau}\int_{0}^{\tau}A(x,y)^{2}(1-G_{c}(x))(1-G_{c^{\prime}}(y))dF_{c}(x)dF_{c^{\prime}}(y)$. Define $J(s_{i},s_{j})$ by

[TABLE]

for all $i,j\in[n]$, and denote by $\tilde{\mathds{E}}$ the expectation conditioned on fixed covariate values $c_{i}$ and $c_{j}$. Then the following hold:

1. i)

$\tilde{\mathds{E}}(J(s_{i},s_{i}))=\int_{0}^{\tau}A(x,x)(1-G_{c_{i}}(x))dF_{c_{i}}(x)$, 2. ii)

$\tilde{\mathds{E}}\left(|J(s_{i},s_{i})|\right)<\infty$, 3. iii)

for all $i\neq j$, $\tilde{\mathds{E}}(J(s_{i},s_{j})|s_{j})=0$, and 4. iv)

for all $i\neq j$, $\tilde{\mathds{E}}\left(J(s_{i},s_{j})^{2}\right)=\int_{0}^{\tau}\int_{0}^{\tau}A(x,y)^{2}(1-G_{c_{i}}(x))(1-G_{c_{j}}(y))dF_{c_{i}}(x)dF_{c_{j}}(y)$.

In our applications of Proposition 23, we will choose

[TABLE]

We continue by proving $n|D_{n}|=o_{p}(1)$. Observe that, by symmetry,

[TABLE]

For the first term in Equation (48), notice that by conditioning on $Q$, we obtain

[TABLE]

where the last equality holds since, without loss of generality, the first $Q$ elements of $\mathcal{D}_{n}$ are chosen from $\mathcal{L}_{0}$ and the $Q$ elements of $\mathcal{D}_{n}^{\prime}$ are chosen from $\mathcal{L}_{1}$. Additionally, notice that by Proposition 23.ii, $\mathds{E}(|J(s_{1},s_{1})||c_{1}=0)<\infty$ and $\mathds{E}(|J(s_{1}^{\prime},s_{1}^{\prime})||c_{1}^{\prime}=1)<\infty$.

Also, note that

[TABLE]

where the last limit holds because $\mathds{V}ar(N_{0})=\eta(1-\eta)n$, and $n_{0}/n\to\eta$. Then, by using the previous result, we deduce

[TABLE]

We continue by proving that the second term in Equation (48) is $o_{p}(1)$. Let $\varepsilon>0$, then

[TABLE]

and

[TABLE]

the first inequality holds because the variables $(s_{i})_{i=1}^{Q}$ are independent of $(s_{i}^{\prime})_{i=1}^{Q}$, given $Q$. The second equality follows from noticing that $\mathbb{C}ov(J(s_{i},s_{j}),J(s_{l},s_{k}))=0$ for all $(i,j)\neq(l,k)$, which is a simple computation that follows straightforwardly from Proposition 23.iii. By Proposition 23.iv, there exists $C^{\prime}>0$ such that for all $i\neq j$, $\mathds{E}(J(s_{i},s_{j})^{2})\leq C^{\prime}$ and $\mathds{E}(J(s_{i}^{\prime},s_{j}^{\prime})^{2})\leq C^{\prime}$, then

[TABLE]

for some $C>0$. We conclude

[TABLE]

where the last limit follows from Equation (49).

Finally, we prove $n|R_{n}|=o_{p}(1)$. Let $\varepsilon>0$, then

[TABLE]

and

[TABLE]

where the previous equality follows by Proposition 23.iii, from which we deduce that all non-trivial covariances are 0, and the inequality holds for some constant $C>0$ due to Proposition 23.iv, as indeed, $\mathds{V}ar\left(J(s_{i},s_{j})|Q\right)$ and $\mathds{V}ar\left(J(s_{i}^{\prime},s_{j}|Q)\right)$ are bounded for any $i,j\in[n]$, independently of the value of the covariates. Since $\mathds{E}(Q)/n\to 0$ we conclude that $nR_{n}=o_{p}(1)$ which, combined with the fact that $nD_{n}=o_{p}(1)$, yields $n(V_{n}-V_{n}^{\prime})=o_{p}(1)$.

∎

B.4 Step III: Limit distribution under the null hypothesis

In this section we prove Theorem 5 for the random covariates model. Recall that by Lemma 22, the result also holds for the deterministic covariates model.

Recall from Equation (45) that

[TABLE]

where $J$ is defined in Equation (44) and $s_{i}=(X_{i},\Delta_{i},c_{i})$ are independent and identically distributed random observations as the covariates are random, which is the main advantage of working in the random covariates model. Then, by using the previous expression, the asymptotic distribution of $(n_{0}n_{1})/nZ_{n}$ can be obtained from a standard application of the theory of V-statistics. In particular, by Proposition 23.iii, $\mathds{E}(J(s_{i},s_{j})|s_{j})=0$ for any $i\neq j$, hence $J$ is a so-called degenerate V-statistic kernel, therefore Theorem 4.3.2 of Koroljuk and Borovskich [1994] yields

[TABLE]

where $\mathcal{Y}=\sum_{i=1}^{\infty}\lambda_{i}(\xi_{i}^{2}-1)$, $(\xi_{i})_{i\geq 1}$ are independent and identically distributed standard normal random variables, and $(\lambda_{i})_{i\geq 1}$ are the eigenvalues associated to the integral operator $T_{J}:L^{2}(X_{1},\Delta_{1},c_{1})\to L^{2}(X_{1},\Delta_{1},c_{1})$ defined by $(T_{J}\phi)(s)=\mathds{E}_{s_{1}}(J(s,s_{1})\phi(s_{1}))$, where $\mathds{E}_{s_{1}}$ means expectation with respect $s_{1}=(X_{1},\Delta_{1},c_{1})$. From Equation (51) we deduce that

[TABLE]

We finish by computing the asymptotic mean and variance of our test statistic. By using Proposition 23.i, a simple computation shows that

[TABLE]

and since $\mathcal{Y}$ has mean 0, the asymptotic mean of $(n_{0}n_{1})/nZ_{n}$ is given by

[TABLE]

Finally, from Proposition 23.iv, we have that

[TABLE]

from which we deduce that the expression for the asymptotic variance is given by

[TABLE]

Appendix C Wild Bootstrap

By following exactly the same steps used in the proof of Theorem 5, the Wild Bootstrap test statistic $Z_{n}^{\mathcal{W}}$, given in Equation (21), can be rewritten as $Z_{n}^{\mathcal{W}}=V_{n}^{\mathcal{W}}+o_{p}(n^{-1})$, where $V_{n}^{\mathcal{W}}=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}W_{i}W_{j}J(s_{i},s_{j})$, and $J$ is the kernel defined in Equation (44). Under the random covariates model $V_{n}$ is a degenerate $V$-statistic, and thus, Theorem 9 is a direct application of Theorem 3.1 of Dehling and Mikosch [1994]. This result can be extended to the deterministic covariates model by using the same coupling of Lemma 22 (however, randomness is taken over $\mathcal{W}$, as we are conditioning on the sequence of data points).

Appendix D Proof of Theorem 6 and Corollary 7

In this section we mainly work under the alternative hypothesis, therefore we need to consider this fact when computing some limit results. For instance, in this more general setting the limit of $n/(n_{0}n_{1})L(x)$ is slightly different to its limit under the null hypothesis. Indeed, by Proposition 10, for every fixed $x$, it holds

[TABLE]

as the number of observations tends to infinity, and by Proposition 11

[TABLE]

uniformly for all $x\leq\tau_{n}$. In the previous results, $\psi^{*}(x)$ is the function defined in Equation (15) which satisfies the following bound

[TABLE]

Other objects that may vary their limits in this more general setting are the pooled Nelson-Aalen estimator $\widehat{\Lambda}$ and Kaplan-Meier statistic $\widehat{F}$, and in general, any statistic constructed using the survival times $(T_{i})_{i\in[n]}$ (as now, $F_{0}$ and $F_{1}$ may or may not be equal). Having this in mind, we proceed to prove Theorem 6 and Corollary 7.

To proof Theorem 6, we require the following intermediate result:

Lemma 24**.**

Assume that Conditions 3 and 4 hold. Then, $\|\phi_{c}^{n}-\phi_{c}\|_{\mathcal{H}}\overset{\mathds{P}}{\to}0$ for $c\in\{0,1\}$.

Proof.

We prove the result for $c=0$, as the proof for $c=1$ is exactly the same. We consider the deterministic covariates model since for the random covariates model we can condition on the number of data points with covariate equal to 0.

Define the intermediate embeddings, $\widehat{\phi}_{0}\in\mathcal{H}$ and $\phi^{\prime}_{0}\in\mathcal{H}$, by

[TABLE]

and recall that $\phi_{0}^{n}$ and $\phi_{0}$ (introduced in Equations (10) and (18), respectively) are given by

[TABLE]

Notice that in this setting the pooled distribution $F$ is not necessarily equal to $F_{0}$.

Clearly

[TABLE]

We will prove that the three terms in the right-hand-side of the previous equation are $o_{p}(1)$. First we prove $\|\phi_{0}^{n}-\widehat{\phi}_{0}\|_{\mathcal{H}}=o_{p}(1)$. Recall that the Nelson-Aalen estimator satisfies $(d\widehat{\Lambda}_{0}-d\Lambda_{0})(x)=dM_{0}(x)/Y_{0}(x)$, then

[TABLE]

where

[TABLE]

We prove that Equation (55) is $o_{p}(1)$ by using Theorem 17. Observe that $h(x,y)$ is symmetric and positive-definite, $(h(x,y))_{(x,y)\in C}$ is $\mathcal{P}$ measurable by Proposition 14, and $(h(x,x))_{x\geq 0}$ is predictable with respect to $(\mathcal{F}_{x})_{x\geq 0}$. Then, by Theorem 17, the result is deduced from checking that $\frac{1}{n_{0}^{2}}\int_{0}^{\tau_{n}}h(x,x)d\langle M_{0}\rangle(x)=o_{p}(1).$ By using that $Y_{1}(x)/Y(x)\leq 1$, $n/n_{1}=O(1)$, and that $Y_{0}(x)/n_{0}=O_{p}(1)S(x)(1-G_{0}(x))$ uniformly for all $x\leq\tau_{n}$, we have

[TABLE]

The previous term is tends to 0 since $\int_{0}^{\tau_{n}}K(\widehat{F}(x-),\widehat{F}(x-))(1-G_{0}(x))dF_{0}(x)<\infty$ by Condition 4.

We continue by proving $\|\widehat{\phi}_{0}-\phi^{\prime}_{0}\|_{\mathcal{H}}=o_{p}(1)$. Let $t<\tau$ and $\epsilon>0$, then

[TABLE]

where $h(x,\cdot)=(K(\widehat{F}(x-),\cdot)-K(F(x),\cdot))$. By taking limsup when $n$ grows to infinity, the first probability is equal to zero since $h(x,y)=\langle h(x,\cdot),h(y,\cdot)\rangle_{\mathcal{H}}\to 0$ uniformly for all $x,y\leq t$. Then,

[TABLE]

Finally, by taking limit when $t\to\tau$, and by using the definition of $h$, we obtain

[TABLE]

We prove that each of the terms in the right-hand-side of the previous Equation are 0. For the first one, from Propositions 11.ii and 11.iii, we have that $(n/(n_{0}n_{1}))L(x)=\beta^{-1}\psi^{*}(x)S_{0}(x)S_{1}(x)$ for all $x\leq\tau_{n}$ with probability at least $1-p(\beta)$ with $p(\beta)\to 0$ as $\beta\to 0$. Then, with probability at least $1-p(\beta)$ we have

[TABLE]

where the second inequality holds by noting that $2K(s,t)\leq K(s,s)+K(t,t)$ for all $s,t$ since $K$ is positive definite, and the fourth inequality holds Equation (54). Finally,

[TABLE]

then, by Condition 4, we conclude that for all $\beta\in(0,1)$

[TABLE]

then taking $\beta\to 0$ yields the desired result. A similar argument shows that

[TABLE]

concluding that $\|\widehat{\phi}_{0}-\phi_{0}^{\prime}\|_{\mathcal{H}}^{2}=o_{p}(1).$

Finally, we prove that $\|\phi_{0}^{\prime}-\phi_{0}\|_{\mathcal{H}}=o_{p}(1)$. Observe that by triangular inequality,

[TABLE]

Observe that

[TABLE]

where $\sigma(x)=\psi^{*}(x)S_{0}(x)S_{1}(x)-n/(n_{0}n_{1})L(x)$. The second equality follows from the dominated convergence theorem: by Equation (52), we have $\sigma(x)\to 0$ for all $x\leq\tau$ almost surely, also, by Equations (53) and (54), it holds $\sigma(x)/S_{0}(x)=O_{p}(1)\psi^{*}(x)S_{1}(x)=O_{p}(1)(1-G_{0}(x))$ uniformly for all $x\leq\tau_{n}$. From Condition 3 we have that $\int_{0}^{\tau}\int_{0}^{\tau}|K(F(x),F(y))|(1-G_{0}(x))(1-G_{0}(y))dF_{0}(x)dF_{0}(y)<\infty$, then we can apply the dominated convergence theorem on sets of probability as high as desired, obtaining that $\|\phi_{0}^{\prime}-\phi_{0}\|_{\mathcal{H}}=o_{p}(1)$.

Finally, observe that

[TABLE]

where the second equality is due to Equation (54), and the third is due Condition 3 together with the fact that $\tau_{n}\to\tau$.

∎

Proof of Theorem 6.

Notice that $\sqrt{Z_{n}}=\|\phi_{0}^{n}-\phi_{1}^{n}\|_{\mathcal{H}}$. Then, by using the triangle inequality, we deduce

[TABLE]

hence, a straightforward application of Lemma 24 deduces $\sqrt{Z_{n}}\overset{\mathds{P}}{\to}\|\phi_{0}-\phi_{1}\|_{\mathcal{H}},$ and by squaring both sides, we obtain the desired result. ∎

Proof of Corollary 7.

Theorem 6 yields $\sqrt{Z_{n}}\overset{\mathds{P}}{\to}\|\phi_{0}-\phi_{1}\|_{\mathcal{H}}$, therefore it is enough to prove that $\phi_{0}\neq\phi_{1}$. Recall from Equation (18) that for $c\in\{0,1\}$, $\phi_{c}$ is defined by

[TABLE]

where the measures $\nu_{c}$ is defined in Equation (17). Let $F^{-1}(x)=\inf\{t\geq 0:F(t)=x\}$, which is differentiable almost everywhere since $F$ is the cumulative distribution function of a continuous random variable. Then, by performing a change of variables, we obtain

[TABLE]

where $f$ is the density function associated with $F$. We conclude that $\phi_{c}$ is the mean kernel embedding (see Equation (6)) of the measure $\pi_{c}(\cdot)$ on $[0,1)$ given by $\pi_{c}(A)=\int_{A}\frac{1}{f(F^{-1}(x))}d\nu_{c}(F^{-1}(x))$.

As the kernel $K$ is continuous and c-universal (recall the definition from Section 3.1), the mean kernel embedding of finite signed measures is injective, hence, we just need to check that $\pi_{0}\neq\pi_{1}$ and that both measures are finite.

First, we prove that $F_{0}\neq F_{1}$ implies that $\pi_{0}\neq\pi_{1}$. We proceed by contradiction, i.e., we assume that $\pi_{0}=\pi_{1}$. Since $F$ and $F^{-1}$ are non-decreasing functions, $\pi_{0}=\pi_{1}$ implies $\nu_{0}=\nu_{1}$. Observe that $\nu_{0}$ and $\nu_{1}$ are continuous measures with respective densities $\psi^{*}(x)S_{1}(x)f_{0}(x)$ and $\psi^{*}(x)S_{0}(x)f_{1}(x)$, where $f_{0}$ and $f_{1}$ denote the density functions of $F_{0}$ and $F_{1}$, respectively. Since $\nu_{0}=\nu_{1}$, their densities are equal. Let $\tau_{0}=\sup\{t:S_{0}(t)>0\}$ and $\tau_{1}=\sup\{t:S_{1}(t)>0\}$, and, without lost of generality, assume that $\tau_{0}\leq\tau_{1}$. By hypothesis, for all $x<\tau_{0}$ we have that $G_{0}(x)<1$ and $G_{1}(x)<1$, hence $\psi^{*}(x)>0$. Using that $\psi^{*}(x)S_{1}(x)f_{0}(x)=\psi^{*}(x)S_{0}(x)f_{1}(x)$, we deduce that for all $x<\tau_{0}$, $S_{1}(x)f_{0}(x)=S_{0}(x)f_{1}(x)$; thus, $\lambda_{0}(x)=\lambda_{1}(x)$, and therefore, by the bijection between hazard and density functions, we have $F_{0}(x)=F_{1}(x)$. Since $F_{0}(\tau_{0})=1$, we conclude that $F_{0}(x)=F_{1}(x)$ for all $x\in\mathds{R}_{+}$, which is a contradiction.

Finally, a simple computation verifies that $\nu_{c}$ is a finite measure. Indeed, for $c=0$, we have $\int_{0}^{\tau}d\nu_{0}(x)=\int_{0}^{\tau}\psi^{*}(x)S_{1}(x)f_{0}(x)dx$, and from Equation (54) it holds that $\psi^{*}(x)S_{1}(x)f_{0}(x)\leq Cf_{0}$, where $C>0$, therefore $\nu_{0}$ is a finite measure. The same holds for $\nu_{1}$, concluding the result.

∎

Appendix E Deferred proofs

E.1 Proof of Theorem 17

To ease notation, we write $h$ and $W$ instead of $h_{n}$ and $W_{n}$, respectively. Define the process $R$ by

[TABLE]

Observe that $R(0)=0$ and that $R(t)\geq 0$ for every $t\geq 0$ since $h$ is positive definite. We will prove that exists an increasing predictable process $A$ that majorises $R$ in the sense that $\mathds{E}(R(T))\leq\mathds{E}(A(T))$ for every finite stopping time, and then, we will apply the Lenglart-Rebolledo inequality to $R$ and $A$, which tells us that for every $\varepsilon,\delta>0$, $\mathds{P}(\sup_{t\leq T}R(t)\geq\varepsilon)\leq\varepsilon/\delta+\mathds{P}(A(T)\geq\delta)$ for any stopping time $T$ (even unbounded).

To find the process $A$, we start by noticing that $R$ can be decomposed as $R(t)=D(t)+2Z(t)$, where the processes $D$ and $Z$ are given by

[TABLE]

respectively. We make two observations. First, since the process $(h(x,y)_{(x,y)\in C}$ is $\mathcal{P}$-measurable, then $(Z(t))_{t\geq 0}$ is an $(\mathcal{F}_{t})$-martingale by Theorem 15. Additionally, by the Optional Stopping Theorem, $Z(T)=0$ for every bounded stopping time $T$. Second, since $h$ is positive definite, $h(t,t)\geq 0$ and $(h(x,x))_{x\geq 0}$ is $(\mathcal{F}_{x})$-predictable, then the process $(D(t))_{t\geq 0}$ is increasing and adapted with $D(0)=0$ and it is compensated by the process $A$ given by

[TABLE]

which is increasing and predictable with respect to $(\mathcal{F}_{t})_{t\geq 0}$. Since $D-A$ is a martingale, we have $\mathds{E}(D(T))=\mathds{E}(A(T))$ for any bounded stopping time $T$. We conclude that for every bounded stopping time $T$, it holds that $\mathds{E}(R(T))=\mathds{E}(A(T))$, i.e., $A$ majorises $R$. Then, if $A(\tau_{n})=o_{p}(1)$, the Lenglart-Rebolledo inequality concludes that $R(\tau_{n})=o_{p}(1)$, and if $A(\tau_{n})=O_{p}(1)$, then the Lenglart-Rebolledo inequality yields $R(\tau_{n})=O_{p}(1)$.

E.2 Proof of Proposition 23

Item i: Observe that

[TABLE]

where $D(t)$ and $Z(t)$ are defined by

[TABLE]

and notice that $d[M^{i}](x)=dN^{i}(x)$ since the survival times are continuous.

We begin by noticing that

[TABLE]

Also, by Theorem 15 and conditioned on the value of the covariate, $Z(t)$ is a zero-mean $(\mathcal{F}_{t})$-martingale. Then the optional stopping time theorem yields $\tilde{\mathds{E}}(Z(X_{i}))=\tilde{\mathds{E}}(Z(0))=0$, and thus

[TABLE]

Item ii: Observe that

[TABLE]

We continue by proving that all the terms on the right-hand-side of the previous equation have positive expectation.

For the first term we have

[TABLE]

where the last inequality follows from the assumptions made in the Proposition. The second term satisfies

[TABLE]

where the last inequality follows from Lemma 1 of Efron and Johnstone [1990] (see also Remark 2 of Akritas [2000]). Finally, the last term satisfies

[TABLE]

where the second equality holds by the symmetry of the functions $A(x,y)$ and $1-H_{c_{i}}(x\vee y)$, and the last inequality holds by Equation (61).

Item iii: Let $i\neq j$ and define the process

[TABLE]

Notice that conditioned on $s_{j}=(X_{j},\Delta_{j},c_{j})$ and $c_{i}$, the process $(Q(t))_{t\geq 0}$ is a zero-mean martingale. Therefore, by the optional stopping theorem we have

[TABLE]

Item iv: Observe that $\tilde{\mathds{E}}(J(s_{i},s_{j})^{2})=\tilde{\mathds{E}}(Q(X_{i})^{2})$ for $i\neq j$, where $(Q(t))_{t\geq 0}$ is the process defined in Equation (62). Observe that, conditioned on $s_{j}=(X_{j},\Delta_{j},c_{j})$ and $c_{i}$, $(Q(t)^{2})_{t\geq 0}$ is an $(\mathcal{F}_{t})$-submartingale with compensator given by

[TABLE]

Then

[TABLE]

and thus

[TABLE]

where the second equality is due to the independence of $s_{i}$ and $s_{j}$ for $i\neq j$, and the third equality follows by noticing that, conditioned on $c_{j}$, $\left(\int_{0}^{t}A(x,y)dM^{j}(y)\right)^{2}$ is an $(\mathcal{F}_{t})$-submartingale for any fixed $x\in\mathds{R}_{+}$ and its compensator is given by $\int_{0}^{t}A(x,y)^{2}Y^{j}(y)d\Lambda_{c_{j}}(y)$.