A kernel- and optimal transport- based test of independence between covariates and right-censored lifetimes

2.1 Right-censored lifetimes

Let T ∈ ℝ ≥ 0 be a lifetime subject to right-censoring, so that we do not observe T directly, but instead observe Z ≔ min { T , C } for some censoring time C ∈ ℝ ≥ 0 , as well as the indicator Δ ≔ 1 { C > T } . We further observe a covariate vector X ∈ ℝ d , where ℝ d is equipped with the Borel sigma algebra. In total, for an i.i.d. sample of size n, we thus obtain the data D ≔ ( ( x i , z i , δ i ) ) i = 1 n ∈ ( ℝ d × ℝ ≥ 0 × { 0 , 1 } ) n .

The main goal of this paper is developing a test of H 0 : X ⫫ T versus H 1 : X ⫫ T based on the sample D. Throughout this paper we will make the following assumptions.

2.2 Independence testing using kernels

Kernel methods have been successfully used for nonparametric independence- and two-sample testing [4], [8]. We now give some of the relevant background in kernel methods.

In [4] it was proposed that the dependence of X and Y could be quantified by the following measure:

Figure 1:

The algorithm for independence testing using HSIC and a permutation test, resulting in a p-value and a rejection decision.

2.3 Optimal transport

Let A = [ a 1 , … , a | A | ] and B = [ b 1 , … , b | B | ] be multisets of length | A | and | B | with a i , b j ∈ ℝ d . Let v A be a distribution vector on A, by which we mean: v A = ( v 1 A , … , v | A | A ) ∈ ℝ | A | so that v i A ≥ 0 and ∑ i = 1 | A | v i A = 1 . Let v B be a distribution vector on B. Let m ( a , b ) be a metric on ℝ d . Then the earth mover’s distance (EMD) between v A and v B is defined as

subject to

Let U be a random variable taking values in A with distribution v A , and U ′ a random variable taking values in B with distribution v B . In the next section we will write ‘let P be the joint distribution that solves the optimal transport problem between U and U ′ ’. By that we mean P ( U = a i , U ′ = b j ) ≔ P i j for i = 1 , … , | A | , j = 1 , … , | B | , where P i j is the matrix that minimizes Eq. (2). Note that P is a valid joint distribution of U and U ′ . Furthermore, for any i with v i A > 0 it holds that P ( U ′ = b j | U = a i ) = P i j / v i A .

3.1 Objective of the algorithm

To use HSIC for an independence test of X and T, one requires the explicit values l ( T i , T j ) for i , j = 1 , … , n . Since for right-censored data T i may be known only in terms of a lower bound, there is no straightforward way to apply these methods to right-censored data. In this section we propose a transformation of the original, censored dataset, into a synthetic dataset consisting of n observed events. The algorithm uses optimal transport and its goal is twofold: first, it should return a dataset to which we can apply a permutation test with test-statistic HSIC, and obtain correct p-values under the null hypothesis; second, it should return a dataset in which the dependence between X and T is similar to the dependence in the original dataset. Indeed, the transformation of the data is of independent interest: we use the standard permutation test with test statistic HSIC/DCOV, but other statistics could be considered. We do believe that the main value of the transformation lies in how it can be straightforwardly combined with permutation testing on the transformed data and we expect it may be difficult to use the transformation for other purposes.

3.2 Definition of the transformation

To define the algorithm, we use the following notation. Let A = [ a 1 , … , a | A | ] be a multiset where a i ∈ ℝ d for i = 1 , … , | A | . We define Uniform ( A ) = ∑ i = 1 | A | δ ( a i ) / | A | , to be the uniform distribution over the elements, where δ ( a ) is the Dirac measure at a. By B = Remove ( a , A ) we set B to be the multiset that remains when one instance of a is removed from A, and by B = Add ( a , A ) we set B to be the multiset that consists of the elements of A, with an added element a. We further assume that we are given a sample D = ( ( x i , z i , δ i ) ) i = 1 n so that z 1 < z 2 < … < z n : that is, we have n distinct event times and they are labeled in increasing order. For a variable a we use a ← b to change the value of a to b. In computing the optimal transport coupling as defined in Section 2.3, we use the Euclidean metric | | x − y | | = ∑ i = 1 d ( x i − y i ) 2 for x ∈ ℝ d . The proposed transformation is defined in Algorithm 2 in Figure 2. Figure 3 visualizes a dataset with 1-dimensional covariates before and after applying the transformation.

Figure 2:

The optHSIC algorithm, resulting in a transformed dataset D ~ .

Figure 3:

Above: The dataset as originally observed. Below: the synthetic dataset after applying the OPT-based transformation. The labels indicate which individual the observation corresponds to in each dataset. The data is sampled from a parabolic relationship between covariates and lifetimes.

3.3 Comments on the transformation algorithm

We now give a verbal explanation of Algorithm 2.

Initialization: The input is simply D = ( ( x i , z i , δ i ) ) i = 1 n where z 1 < … < z n . We initialize an empty transformed dataset D ˜ , to which we will add n observations of the form ( x i , t j ) for some i , j ∈ { 1 , … , n } in the following way. We will loop over the times z i from i = 1 , … , n . At each time, the multiset AR lists the covariates at risk in the dataset D and AR is initialized as the multiset containing all covariates. The multiset L will list all covariates that have not been added to D ˜ yet. Indeed, as D ˜ is initialized empty, L initially contains all n covariates. The variable d will count the number of observed δ = 1 events and is initialized 0.

First loop from i = 1 , … , n − 1 : At time z i we distinguish two cases: if δ i = 0 we leave L and D ˜ unchanged, and simply remove one instance of x i from AR. If δ i = 1 we add 1 to d (to count the number of observed events). We also select a covariate x ˜ from L as follows: First a joint distribution of Uniform ( AR ) and Uniform ( L ) is computed using optimal transport. This is a matrix P of size | AR | × | L | . This distribution is then conditioned on the event that Uniform ( AR ) = x i , yielding a distribution over L. We sample from this distribution to obtain the covariate x ˜ . The pair ( x ˜ , z i ) is then added to D ˜ , and the covariate x ˜ is removed from L. There are now d observations in D ˜ and there are n − d covariates left in L. We also remove x i from AR.

When we finish the first loop: It is now the case that AR = [ x n ] (note the previous loop went up to i = n − 1 ) and the multiset L contains n − d covariates. In the second loop, that runs from d + 1 to n, each remaining covariate in L is combined with the time z n and added to D ˜ . After this loop, L is empty, and D ˜ contains all n covariates x 1 , … , x n , each associated with a time. The transformed dataset D ˜ is returned.

Consider the special case in which there is no censoring and δ i = 1 for i = 1 , … , n . Then it is easy to verify that in the first loop AR = L at each step i, so that the optimal transport algorithm returns a scaled identity matrix. Consequently, at step i the covariate x ˜ is equal to x i with probability 1, and the final output is D ˜ = D . A nontrivial example is worked out in Section A.1.

3.4 Intuition behind the transformation

Before we prove properties of the proposed transformation, we briefly comment on the intuition behind the transformation. To this end, first consider a permutation test in the absence of censoring, when we simply observe D = ( ( x i , t i ) ) i = 1 n . The permuted datasets can be generated as follows: loop through the events in order of time, and to each time, associate a covariate that you have not associated to any earlier time. Comparing the original dataset with the datasets obtained in this manner is justified for the following reason: Under the null hypothesis a sample can be generated by firstly sampling t i for i = 1 , … n i.i.d. and independently sampling x i for i = 1 , … n i.i.d., and secondly looping through the times in order, and associating to each time a covariate that has not yet been chosen at a previous time. The original dataset and the permuted datasets are thus equal in distribution: intuitively, the permutation test checks whether the dataset looks as if, at each time, a covariate is picked uniformly from those not chosen before.

It is not obvious how to translate this to censored data. Due to censoring, it may not be true that the ith event covariate is chosen uniformly from the covariates that have not had an observed event time before. For this reason survival analysis often compares the i–th event covariate with the covariates at risk (not failed and not censored) just before the i–th event. If the null hypothesis holds, then intuitively it holds that the i–th event covariate is chosen uniformly from the covariates at risk at time z i . That is, if AR denotes the multiset of covariates at risk and x i is the event covariate, we would like to test if x i was chosen uniformly from AR.

To do so using a permutation test, our algorithm couples Uniform ( AR ) to a uniform choice from those that have not yet been chosen in the synthetic dataset: Uniform ( L ) (see Algorithm 2). In this manner, when the null hypothesis is true, (intuitively) it holds that in the synthetic dataset the covariates are chosen uniformly from those not chosen before. But this last statement is exactly our intuition behind a permutation test: namely, we use a permutation test to see if the dataset looks as if, at each time, a covariate is picked uniformly from those not chosen before.

The intuition discussed thus far related to the workings of our transformation under the null hypothesis, and had at its core that Uniform ( AR ) is coupled to Uniform ( L ) in Algorithm 2. However there are many couplings between these two distributions. The reason we choose the coupling that solves the optimal transport problem is the following: assume the alternative hypothesis holds and that given AR _i , certain covariates have a higher hazard rate than others. Then we would like this bias towards these covariates to be visible in D ˜ . In other words we would like x ˜ to be close to x ˜ i in Algorithm 2. That is precisely what the optimal transport coupling achieves.

Finally, including the remaining covariates in the transformed dataset at the last time ensures the permuted datasets correctly reflect all alternative choices that can be made when covariates are chosen at random from those not chosen before. The association of these remaining covariates to the last event time reflects they have not been selected at each of the earlier times, which may indicate they have had less risk of having an event up to that point.

The goal of the transformation is thus twofold: firstly, ensuring that under the null hypothesis a permutation test is appropriate on the transformed dataset, which motivates the coupling of Uniform ( AR ) and Uniform ( L ) ; secondly, given the first goal, ensure covariates in D ˜ associated to times z i for δ i = 1 are as close as possible to the covariates associated to z i in D, which motivates the chosen coupling to be the coupling that solves the optimal transport problem defined in Section 2.3.

An interesting alternative approach would be to compare x i with Uniform ( AR i ) directly, without first transforming the data. We do not see, however, how to combine that approach with a permutation test, or with another procedure to test for significance.

4.1 Computational cost of optHSIC

The algorithm optHSIC consists of two parts: First, the dataset is transformed. Second, a permutation test with test statistic HSIC/DCOV is performed. The computational cost of the transformation may be high: the solution of the earth mover’s distance has n 3 log ( n ) complexity [13]. Since the earth mover’s distance is computed for each i = 1 , … , n such that δ i = 1 , complexity of the transformation is excessively high. Luckily, fast approximations of the earth mover’s distance exist that are linear in time (e.g. [13]), making the transformation O ( n 2 ) . Also in the case where the covariates are 1-dimensional, the optimal transport problem has a simple solution [14]. Different algorithms and/or approximations also exist for the EMD with other metrics [13]. The second step, computation of HSIC is also an O ( n 2 ) -operation for which large-scale approximations can be made [15]. Furthermore, both the computation of HSIC and the permutation test computations can be easily parallelized. In our simulations we did not use approximations of the earth mover’s distance nor of HSIC, as for samples up to about 500 the optHSIC test can be performed in about a second on an ordinary PC.

4.2 Theoretical results on optHSIC

We say a test with p-value p has correct type 1 error rate if under the null hypothesis P H 0 ( p ≤ α ) ≤ α . The main theoretical result we obtained for optHSIC is that the type 1 error rate is correct when C ⫫ X . Unfortunately, we have not been able to prove other important results such as (asymptotically) correct type 1 error rate when C ⫫ X , or consistency (power converging to one for each alternative hypothesis). However, as Section 6 details, extensive simulations demonstrate that optHSIC achieves correct type 1 error rate also when C ⫫ X . Simulations further show that optHSIC is able to detect a wide range of dependencies between X and T without losing much power, relative to the CPH likelihood ratio test, even when the CPH model assumptions holds. Obtaining further theoretical results is an important future challenge. In the uncensored case a permutation test with covariance kernel of Brownian motion is consistent for distributions with compact support and has correct type 1 error rate [12]. The difficulty of extending these proofs to optHSIC is that optHSIC requires sequential analysis of the optimal transport distributions, breaking independence between observations in the transformed dataset.

We first prove an auxiliary result: namely, although we propose to permute the transformed dataset, this is equivalent to permuting the original dataset, and then transforming the permuted datasets.

4.3 Using multiple transformations

Algorithm 2 defines the transformation used in optHSIC. In line 7, a covariate is sampled from the multiset L with distribution v. The sampled element may differ for different iterations of the algorithm. Consequently, given a fixed dataset D, the transformed dataset will look different for different iterations of the algorithm. The optHSIC algorithm proposes to use only one transformation, resulting in a single p-value.

Say that instead of applying the optHSIC algorithm once, we repeat the optHSIC algorithm m times, resulting in m p-values p 1 , … , p m . An example of the distribution of p-values for a fixed dataset is given in Figure 5.

Figure 5:

A histogram of the p-values obtained from optHSIC for a fixed dataset D sampled from D.1 of Table 2. The sample was of size n = 300 and 61 % of the individuals was observed ( δ = 1 ) . The variability in p comes from the fact that the transformation in optHSIC is random.

Figure 5 shows that in the used dataset, where 39 % of observations is censored and n = 300 , the variance in the obtained p-values is not excessively high, and in particular all p-values would have led to the decision to reject the null-hypothesis.

We also studied several techniques one may use to combine p-values, typically at the cost of being more conservative when there is little variance among p-values, like in the example above. This is discussed in Section A.9.

5.1 wHSIC

HSIC relies on estimating the mean embedding of the joint distribution P X T . The estimated embedding is then compared to the embedding of the product of marginal distributions. In the uncensored case the empirical distribution equals ∑ i = 1 n δ ( x i , t i ) / n which has corresponding mean embedding ∑ i = 1 n K ( ( x i , t i ) ) / n .

Since we do not observe the t _i ’s we could consider replacing the empirical distribution by a weighted version ∑ i = 1 n w i δ ( x i , z i ) / n where we try to find weights w i such that ∑ i = 1 n w i f ( x i , z i ) ≈ E f ( X , T ) for every measurable function f such that the expectation exists. A natural idea is to give an observed point ( x , z , δ ) a weight of zero if it is censored, and a weight equal to the inverse probability of being uncensored otherwise. This can be motivated by the following lemma, that applies also if C ⫫ X . We first define

The following lemma proposes a weight function W in terms of g.

Figure 6:

The wHSIC algorithm, resulting in a p-value and rejection decision of H 0 : X ⫫ T .

5.2 zHSIC

If C ⫫ X , then any dependence between X and Z must be due to dependence between X and T. Hence we may estimate HSIC ( X , Z ) to measure the strength of the dependence. This motivates the test proposed in Algorithm 5 labeled zHSIC (see Figure 7). This test is, indeed, expected to yield false rejections if C ⫫ X fails to hold, and to lack power when a large portion of the events is censored. We include this test mainly to see how the power of optHSIC and wHSIC compare with this approach.

Figure 7:

The zHSIC algorithm, resulting in a p-value and rejection decision of H 0 : X ⫫ T .

6.1 Type 1 error rate

We begin by investigating the type 1 error rate. The distributions in which we test the type 1 error rate are found in Table 3. In these distributions the null hypothesis holds, i.e. X ⫫ T and we consider both cases where C ⫫ X and C ⫫ X , as well as the case of multidimensional covariates. We estimate the rejection rates by sampling 5000 times from each distribution for each sample size and applying the different tests to the samples. The obtained rejection rates of optHSIC are found in Table 4. The type 1 error rate of the remaining methods is found in Section A.10.1. The most important finding is that optHSIC is found to have correct type 1 error rate of α = 0.05 both when C ⫫ X and when C ⫫ X . In contrast, wHSIC and zHSIC yield many false rejections when C ⫫ X , as expected, but have the correct type 1 error rate when C ⫫ X . Investigation of the p-values of optHSIC furthermore showed that p-values are distributed approximately according to Uniform [ 1 B + 1 , … , B + 1 B + 1 ] under the null hypothesis, even when C ⫫ X (see Figure 11 in Section A.10.1).

6.2 Comparison of power

To compare the power of the four methods we consider a number of distributions in which T ⫫ X . For each distribution we let the sample size range from n = 40 to n = 400 in steps of 40. At each sample size we take 1000 samples to estimate the rejection rate. In these distributions censoring is independent of the covariate such that the methods wHSIC and zHSIC do not have inflated rejection rate due to dependencies between C and X. Distributions with varying censoring rates and dependent censoring are investigated in Section A.10.2. Parameters are chosen so that that 60 % of the observations is uncensored ( 60 % has δ i = 1 ), and rejection rates take a range of values (i.e. to exclude trivial distributions so that each method rejects with probability one for each sample size).

6.2.1 Power for one-dimensional covariates

In distributions 1–4 of Table 2 the covariate is one-dimensional. Scatterplots of the samples and rejection rates are displayed in Figure 8. Note how in D.1 the CPH assumption holds, so the CPH method suits this example very well. We find however that the rejection rate of optHSIC is very similar to that of the CPH likelihood ratio test (first row, right of Figure 8). D.2 features a case in which hazard is highest in the middle, and lower for extreme values of the covariate. The CPH likelihood ratio test does not have power to detect this relationship, and optHSIC is the most powerful method. In D.3 and D.4 the hazard functions of different covariates cross each other. In this case, wHSIC is the top-performing method, but optHSIC is also able to detect the more complicated relationship, while CPH is not able to do so.

Figure 8:

Scatterplots and rejection rates for D.1 (top row) – D.4 (bottom row) of Table 2.

6.2.2 Power for multidimensional covariates

In D.5–8 of Table 2 the covariates are multidimensional. Figure 9 shows the rejection rates of the four methods, both as the dimension is fixed at 10, and the sample size increases (left column) and when the sample size is fixed at 200, and the dimension increases (right column). In D.5 and D.6 the CPH assumption holds. Again we see that the power of optHSIC is relatively close to the power of the CPH likelihood ratio test, which is specifically designed for these assumptions. This holds both when the dependence is on a single sub-dimension of the covariate (D.6) and when the dependence is on all covariates (D.5). In D.7 and D.8 a non-linear term is present in the hazard rate. Together with zHSIC, optHSIC is now the best performing method. Note that as dimension increases, the dependence in D.6–D.8 becomes harder to detect, whereas the dependence in D.5 becomes easier to detect since in the latter T depends on all covariates.

Figure 9:

Rejection rates of the four methods for D.5 (top row) – D.8 (bottom row) of Table 2. Left column: dimension fixed at 10, and the sample size increases. Right column: sample size fixed at 200, and the dimension increases.

We also investigated the power when C ⫫ X and as the censoring percentage varies from 0 % to 80 % . Distributions and rejection rates are presented in Section A.10.2. The main observations are that, although all methods lose power as the censoring rate increases, the relative performance of the four methods remains similar to the results presented in the main text.

We thus find that for continuous covariates optHSIC is able to detect a wider range of dependencies than the CPH likelihood ratio test, while not losing much power when the CPH assumptions hold. This is true both when the covariate is 1-dimensional and when the covariate is multidimensional, even when the dependence arises from a lower-dimensional subspace of the covariate. Importantly, unlike the methods wHSIC and zHSIC, the type 1 error rate of optHSIC is found to be of the correct level both when C ⫫ X and when C ⫫ X .

6.3 Binary covariates

Lastly, we did experiments with binary covariates, in which case independence testing is equivalent to two-sample testing. For two-sample testing there are other alternative approaches, including [17], [18]. Furthermore, the wHSIC approach can be modified firstly to allow for the computation of weights within groups, and secondly a permutation test can be constructed that is valid also when censoring differs between groups. This permutation test was proposed by [19]. This special case of wHSIC coincides with the approach discussed in [18] for two-sample testing.

Our main finding is that optHSIC has a decent all-round performance for two-sample testing. However, especially when the distributions meet additional assumptions used in semi-parametric tests, these semi-parametric tests outperform optHSIC. We also note that, as optHSIC relies on choosing ‘similar covariates’ during the transformation of Algorithm 2, it is not ideally suited for the two-sample case, as there is a chance of choosing the opposite covariate.

We also provide a ‘failure mode’ of optHSIC (we thank a reviewer for this example): a scenario in which optHSIC performs much worse than the CPH test. In this example, covariates are binary (i.e. there are two groups), group sizes are very unequal, 90 % of observations is censored, censoring appears only in the large group, and survival in the two groups is identical until beyond the censoring has occurred. See A.12 for a detailed description.

A full presentation of the experiments and methods used for the case of binary covariates is given in Section A.11. We conclude that the main value of optHSIC lies in independence testing with continuous covariates.