K -sample omnibus non-proportional hazards tests based on right-censored data

Abstract

This work presents novel and powerful tests for comparing non-proportional hazard functions, based on sample–space partitions. Right censoring introduces two major difficulties, which make the existing sample–space partition tests for uncensored data non-applicable: (i) the actual event times of censored observations are unknown and (ii) the standard permutation procedure is invalid in case the censoring distributions of the groups are unequal. We overcome these two obstacles, introduce invariant tests, and prove their consistency. Extensive simulations reveal that under non-proportional alternatives, the proposed tests are often of higher power compared with existing popular tests for non-proportional hazards. Efficient implementation of our tests is available in the R package KONPsurv, which can be freely downloaded from CRAN.

Keywords

Combined test consistent test crossing hazards non-parametric test permutation test right censoring sample–space partition versatile test

1 Introduction

For the task of comparing survival distributions of two or more groups using censored data, the logrank test is the most popular choice. Its optimality properties under proportional-hazard functions are well known. Although the logrank test is asymptotically valid, it may not be powerful when the proportional hazards assumption does not hold. There are varieties of situations in which the hazard functions are of non-proportional shape. For example, a medical treatment might have adverse effects in the short run, yet effective in the long run, or a treatment may be short-term beneficial but gradually lose its effect with time. In such scenarios, the hazard functions cross. In general, the longer the follow-up period is, the more likely it is for various non-proportional scenarios to develop.¹

Other tests have been proposed that might be better choices for non-proportional hazards under the alternative. Peto and Peto² proposed a test, which is similar to the logrank test, but more sensitive for differences in hazards at early survival times than at late ones. Pepe and Fleming^3,4 suggested a weighted Kaplan–Meier (KM) test with a weight function consists of the geometric average of the two censoring survival-function estimators. Yang and Prentice¹ recently proposed another weighted logrank test whose weights are obtained by fitting their model,⁵ which includes the proportional hazards and the proportional odds models as special cases. In contrast to the logrank and other related tests, the test of Yang and Prentice¹ uses adaptive weights. Under proportional hazards alternatives, this new adaptively weighted logrank test is optimal. When the hazards are non-proportional, the adaptive weights typically lead to improvement in power over the logrank test. The test of Yang and Prentice,¹ referred here as the Yang–Prentice test, is currently considered to be the leading one in terms of power, under a wide range of non-proportional hazards alternatives. However, this test is applicable only for two-sample problems. Moreover, it is not invariant to group labeling. Exchanging the group labels between treatment and control would result in a different p-value. Thus, in applications with no clear link between the groups to treatment/control labeling, such as in testing the differences between females and males, the Yang–Prentice test in its current form is inappropriate. In Section 3.2., we suggest an invariant version of the Yang–Prentice test.

In the statistical literature of K-sample tests for non-censored data, there exist powerful consistent tests that are based on various sample–space partitions. These include the well-known Kolmogorv–Smirnov and Cramer–von Mises tests,⁶ and the Anderson–Darling (AD) family of statistics.^7,8 In particular, Thas and Ottoy⁹ showed that the K-sample AD test is basically an average of Pearson statistics in $2 \times K$ contingency tables that are induced by observation-based partitions of the sample space into two subsets. They suggested an extension of the AD test, by considering higher partitions, up to four. Heller et al.¹⁰ proposed the HHG test, a sample–space partition-based non-parametric test for detecting associations between two random vectors of any dimension. When one of the random vectors is a categorical one-dimensional variable, the problem reduces to the K-sample problem with an observation-based partition of the sample space into two subsets, using three intervals. This specific partition is adopted in this work and will be described in detail. Heller et al.¹¹ extended the work of Thas and Ottoy⁹ by considering sample–space partitions higher than 4, and studied test statistics that aggregate over all partitions by summation or maximization and also by aggregating over different sizes of partitions. They showed, by extensive simulation studies, that increasing the number of partitions can increase power under complex settings in which the density functions intersect four times or more.

In this work, we present new powerful non-parametric and invariant tests for comparing two or more survival distributions using right-censored data. Our proposed methodology is demonstrated and applied using the specific sample–space partition of Heller et al.,¹⁰ which has been shown to be very powerful with 3 or less densities’ intersections,¹¹ under non-censored data. Right-censored data introduce two major difficulties: (i) the actual event times of censored observations are unknown; and (ii) the standard permutation procedure of label shuffling is invalid, in case the censoring distributions of the groups are unequal. We overcome these two obstacles and introduce novel consistent powerful tests for right-censored data. Additionally, we provide a robust powerful test under non-proportional or proportional hazards, based on the principle of minimum p-value and the Cauchy-combination test of Liu and Xie.¹² The power of our robust test is compared with the test of Lee,¹³ which is based on the maximum of two weighted logrank test statistics, and the MaxCombo test¹⁴ based on the maximum of the logrank and three weighted logrank test statistics.

2 K-sample tests based on sample–space partition

2.1 Motivation and notation

Let X be a one-dimensional non-negative random variable, $X \in R^{+}$ , and let Y be a categorical variable indicating the group label. Under the K-sample hypotheses testing, the null hypothesis is $H_{0} : F_{1} (x) = \dots = F_{K} (x)$ for all $x \in R^{+}$ , and the alternative is $H_{1} : F_{m} (x) \neq F_{k} (x)$ for some $1 \leq m < k \leq K$ and some $x \in R^{+}$ , where F_k is the cumulative distribution function of group k, $k = 1, \dots, K$ . We assume that the sample spaces, on which these K distributions are defined, coincide.

K random samples $A_{1}, \dots, A_{K}$ are drawn from the respective distributions $F_{1}, \dots, F_{K}$ . Let n_k be the total number of samples in group k, $k = 1, \dots, K$ , and $n = \sum_{k = 1}^{K} n_{k}$ . Assume temporarily no censoring, and consider a pair of observations $X_{i} \in A_{k}$ and X_j. Then, in the spirit of the HHG test,¹⁰ we consider the following partition induced by the pair (i, j): $A_{11} (i, j)$ is the number of observations from group k that the distance between the value of X and X_i is less than or equal to $| X_{i} - X_{j} |$ ; $A_{12} (i, j)$ is the number of observations outside group k that the distance to X_i is less than or equal to $| X_{i} - X_{j} |$ ; $A_{21} (i, j)$ is the number of observations from group k that the distance to X_i is larger than $| X_{i} - X_{j} |$ ; and $A_{22} (i, j)$ is the number of observations outside group k that the distance to X_i is larger than $| X_{i} - X_{j} |$ .

Interestingly, with no censoring, the A_lr’s, $l, r = 1, 2$ , can be expressed through the empirical cumulative distribution functions. For example, assume X_i > X_j and $Y_{i} = Y_{j} = k$ . Then

\begin{array}{l} A_{11} (i, j) = \sum_{r = 1, r \neq i, j}^{n} I {| X_{i} - X_{r} | \leq X_{i} - X_{j}} I {Y_{r} = k} \\ = \sum_{r = 1, r \neq i, j}^{n} I {X_{j} \leq X_{r} \leq 2 X_{i} - X_{j}} I {Y_{r} = k} \\ = n_{k} {{\hat{F}}_{k} (2 X_{i} - X_{j}) - {\hat{F}}_{k} (X_{j}^{-})} - 2 \end{array}

where

{\hat{F}}_{k} (x) = n_{k}^{- 1} \sum_{X_{l} \in A_{k}} I (X_{l} \leq x)

and

{\hat{F}}_{k} (x^{-}) = n_{k}^{- 1} \sum_{X_{l} \in A_{k}} I (X_{l} < x)

. The –2 above stands for excluding observations i and j. In general, all observations with a distances to X_i that is less than or equal to

| X_{i} - X_{j} |

lie inside the interval

[a_{i j}, b_{i j}]

, where

a_{i j} = \min (X_{j}, 2 X_{i} - X_{j})

and

b_{i j} = \max (X_{j}, 2 X_{i} - X_{j})

, as illustrated in Figure 1.

Figure 1.

sample–space partition.

In general, for Y _i = k, we get

\begin{array}{l} A_{11} (i, j) = n_{k} {{\hat{F}}_{k} (b_{i j}) - {\hat{F}}_{k} (a_{i j}^{-})} - 1 - I (Y_{i} = Y_{j}), \\ A_{12} (i, j) = \sum_{m = 1, m \neq k}^{K} n_{m} {{\hat{F}}_{m} (b_{i j}) - {\hat{F}}_{m} (a_{i j}^{-})} - I (Y_{i} \neq Y_{j}), \\ A_{21} (i, j) = n_{k} - A_{11} (i, j) - 1 - I (Y_{i} = Y_{j}), \\ A_{22} (i, j) = \sum_{m = 1, m \neq k}^{K} n_{m} - A_{12} (i, j) - I (Y_{i} \neq Y_{j}) \end{array}

For each pair (i, j), a 2 × 2 contingency table can be constructed with $A_{l r} (i, j)$ as the entry of cell lr, $l, r = 1, 2$ , and a total sum of n – 2. Under the null hypothesis of equal distributions, the probability of belonging to cell lr equals the product of the marginal probabilities of the lth row and the rth column. Therefore, a summary statistic of a such contingency table can be based on either the Pearson chi-squared test statistic, or the log-likelihood ratio statistic. Since it is unknown which pair (i, j) yields the best sample–space partition that provides the largest summary statistic, the final test statistic is defined as a sum over all possible $n (n - 1)$ partitions induced by the data. A permutation-based p-value can be calculated, under random permutations of the group labels.

To introduce the right-censored data, let $C \in R^{+}$ be a non-negative random variable indicating the censoring time. Assume that X and C are conditionally independent given Y. Define T to be the observed time, namely $T = \min (X, C)$ and let $Δ = I (X \leq C)$ . Hence, the observed data consist of K random samples that can be summarized by $(T_{i}, Δ_{i}, Y_{i}), i = 1, \dots, n$ . Note that the different groups may have different censoring distributions. With right-censored data, our proposed test requires special care in evaluating the $A_{r, l}$ 's, $r, l = 1, 2$ , by replacing the empirical distribution functions by their respective KM estimators, as well as when applying a permutation test with unequal censoring distributions. Both issues are described in details below.

2.2 The test statistic

Let ${\tilde{F}}_{k}$ be the KM estimator of the cumulative distribution function using all observations of group k. The KM estimator is defined only up to (and including) the last observed failure time. Define γ_k to be the maximum time in which ${\tilde{F}}_{k}$ can be used for the test statistic. If the largest observed time is a failure time, ${\tilde{F}}_{k} (t)$ is known for the entire range of t. In this case, we define γ_k to be the maximum possible value of t required for the test, $2 \max_{i = 1, \dots, n} {T_{i}} - \min_{i = 1, \dots, n} {T_{i}}$ . However, in case of censoring after the largest failure time, the KM estimator beyond the largest observed failure time is undefined. We thus define the maximum time γ_k to be the largest observed failure time. Namely

γ_{k} = {\begin{array}{l} 2 \max_{i = 1, \dots, n} {T_{i}} - \min_{i = 1, \dots, n} {T_{i}} if \max_{i = 1, \dots, n} {T_{i} Δ_{i} I (Y_{i} = k)} = \max_{i = 1, \dots, n} {T_{i} I (Y_{i} = k)} \\ \max_{i = 1, \dots, n} {T_{i} Δ_{i} I (Y_{i} = k)} if \max_{i = 1, \dots, n} {T_{i} Δ_{i} I (Y_{i} = k)} \neq \max_{i = 1, \dots, n} {T_{i} I (Y_{i} = k)} \end{array}

Define $γ_{- k}$ to be the maximum time in which at least one of the other KM estimators, ${\tilde{F}}_{m}, m = 1, \dots, K, m \neq k$ , can be used for the test statistic, namely $γ_{- k} = \max_{m \neq k} {γ_{m}}$ . Define τ_k to be the maximum time point the KM estimator is defined in group k and in at least one more group, namely $τ_{k} = \min {γ_{k}, γ_{- k}}$ . Then, for each pair of observed failure times $T_{i} \in A_{k}$ and T_j such that $j \neq i, Δ_{i} = Δ_{j} = 1$ and $b_{i j} \leq τ_{k}$ , a 2 × 2 contingency table is constructed. The following $A_{l r}^{*} (i, j), l, r = 1, 2$ , are the corresponding expressions of $A_{l r} (i, j)$ obtained by replacing $\hat{F}$ by $\tilde{F}$ , the KM estimators. In case $γ_{m} < b_{i j}$ , the observations of group m will not be included in the contingency table induced by the pair (i, j). Namely, for $j \neq i, Δ_{i} = Δ_{j} = 1, b_{i j} \leq τ_{k}$ and Y _i = k

\begin{array}{l} A_{11}^{*} (i, j) = n_{k} {{\tilde{F}}_{k} (b_{i j}) - {\tilde{F}}_{k} (a_{i j}^{-})} - 1 - I (Y_{i} = Y_{j}), \\ A_{12}^{*} (i, j) = \sum_{m = 1, m \neq k}^{K} n_{m} {{\tilde{F}}_{m} (b_{i j}) - {\tilde{F}}_{m} (a_{i j}^{-})} I (γ_{m} \geq b_{i j}) - I (Y_{i} \neq Y_{j}), \\ A_{21}^{*} (i, j) = n_{k} - A_{11}^{*} (i, j) - 1 - I (Y_{i} = Y_{j}), \\ A_{22}^{*} (i, j) = \sum_{m = 1, m \neq k}^{K} n_{m} I (γ_{m} \geq b_{i j}) - A_{12}^{*} (i, j) - I (Y_{i} \neq Y_{j}) \end{array}

Only pairs of observed failure times are used for the sample–space partitioning (i.e. $Δ_{i} = Δ_{j} = 1$ ), while censored observations contribute in ${\tilde{F}}_{k}, k = 1, \dots, K$ . Denote by n(i, j) the number of observations in all the groups included in the contingency table induced by (i, j), namely $n (i, j) = \sum_{k = 1}^{K} n_{k} I (γ_{k} \geq b_{i j})$ . Let S_P and S_LR be the summary statistic of each contingency table, based on Pearson chi-squared test statistic

S_{P} (i, j) = \frac{{n (i, j) - 2} {A_{12}^{*} (i, j) A_{21}^{*} (i, j) - A_{11}^{*} (i, j) A_{22}^{*} (i, j)}^{2}}{A_{1 \cdot}^{*} (i, j) A_{2 \cdot}^{*} (i, j) A_{\cdot 1}^{*} (i, j) A_{\cdot 2}^{*} (i, j)}

and the log-likelihood ratio statistic

S_{L R} (i, j) = 2 \sum_{m = 1, 2} \sum_{k = 1, 2} A_{m k}^{*} (i, j) \log [\frac{{n (i, j) - 2} A_{k l}^{*} (i, j)}{A_{m \cdot}^{*} (i, j) A_{\cdot k}^{*} (i, j)}]

respectively, where

A_{. k}^{*} (i, j) = \sum_{m = 1, 2} A_{m k}^{*} (i, j)

and

A_{m .}^{*} (i, j) = \sum_{k = 1, 2} A_{m k}^{*} (i, j)

. In case of at least one zero margin in the contingency table,

S_{P} (i, j) = 0

and

S_{L R} (i, j) = 0

. Denote by

N = \sum_{k = 1}^{K} \sum_{i = 1, T_{i} \in A_{k}}^{n} \sum_{j = 1, j \neq i}^{n} Δ_{i} Δ_{j} I (T_{j} \leq τ_{k}) I (b_{i j} \leq τ_{k})

the total number of tables constructed from the data. Then, our proposed sample–space partition test statistic for equality of K distributions based on right-censored data is defined by

Q = \frac{1}{N} \sum_{k = 1}^{K} \sum_{i = 1, T_{i} \in A_{k}}^{n} \sum_{j = 1, j \neq i}^{n} S (i, j) Δ_{i} Δ_{j} I (b_{i j} \leq τ_{k})

where S(i, j) is either the test statistic

S_{P} (i, j)

S_{L R} (i, j)

. In the case of no right censoring, the number of tables is solely determined by the number of observations. In right-censored data, the number of tables is random and determined also by the actual observed values due to the restrictions

Δ_{i} = Δ_{j} = 1

and

b_{i j} \leq τ_{k}

. This issue is of high importance for the permutation stage of the test, as elaborated in the next subsection.

2.3 The permutation procedure

Allegedly, a permuted test can be done based on random permutations of the group labels. However, if the censoring distributions of the K groups are different, such a permutation test is invalid, since a significant result can be yielded under the null due to differences in the censoring distributions. In order to generate random permutations that are independent of Y, we adopt the imputation approach suggested by Wang et al.¹⁵

The main idea consists of randomly permuting the group labels, while for each observation assigned to a group different from the original one; a censoring time is imputed from the censoring distribution of the new assigned group. If the observation was originally censored, a survival time is also imputed, from the null survival distribution. Let $Y_{1}^{p}, \dots, Y_{n}^{p}$ be a random permutation of the group labels. Define $(T_{i}^{p}, Δ_{i}^{p}), i = 1, \dots, n,$ by

(T_{i}^{p}, Δ_{i}^{p}) = {\begin{array}{l} (T_{i}, Δ_{i}) if Y_{i}^{p} = Y_{i} \\ ({\tilde{T}}_{i}, {\tilde{Δ}}_{i}) if Y_{i}^{p} \neq Y_{i} \end{array}

where

{\tilde{T}}_{i} = \min ({\tilde{X}}_{i}, {\tilde{C}}_{i}), {\tilde{Δ}}_{i} = I ({\tilde{X}}_{i} \leq {\tilde{C}}_{i})

. Furthermore,

{\tilde{C}}_{i}

is sampled from the estimated censoring distribution of group

Y_{i}^{p}

, based on the KM estimator of the censoring distribution of group

Y_{i}^{p}

, by reversing the roles of event and censoring. The KM estimator, denoted by

{\hat{G}}_{C, Y_{i}^{p}}

, is defined up to the largest observed censoring value of that group. Then, each observed censoring time is sampled with probability equals to the jump size of the respective KM estimator. In case

{\hat{G}}_{C, Y_{i}^{p}}

is incomplete, i.e.

{\hat{G}}_{C, Y_{i}^{p}} [\max_{j = 1, \dots, n} {T_{j} (1 - Δ_{j}) I (Y_{j} = Y_{i}^{p})}] > 0

, the largest observed time,

\max_{j = 1, \dots, n} {T_{j} (1 - Δ_{j}) I (Y_{j} = Y_{i}^{p})}

, is sampled with probability

{\hat{G}}_{C, Y_{i}^{p}} [\max_{j = 1, \dots, n} {T_{j} (1 - Δ_{j}) I (Y_{j} = Y_{i}^{p})}]

{\tilde{X}}_{i}

is defined by

{\tilde{X}}_{i} = {\begin{array}{l} X_{i} if Δ_{i} = 1 \\ X_{i}^{*} if Δ_{i} = 0 \end{array}

where

X_{i}^{*}

is sampled from an estimator of

pr (X_{i} > x | X_{i} > T_{i})

, the conditional distribution of X under the null hypothesis. In practice,

pr (X_{i} > x | X_{i} > T_{i})

is replaced by its KM estimator, using all observations from all groups that their observed time is larger than T_i. Denote this KM estimator by

{\hat{S}}_{cond, T_{i}}

. The sampling based on

{\hat{S}}_{cond, T_{i}}

is done in the same fashion as in the above censoring sampling, but in case of incomplete distribution, the value

\max_{i = 1, \dots, n} (T_{i} Δ_{i}) + ε

(ε is any positive number) is sampled with probability

{\hat{S}}_{cond, T_{i}} {\max_{i = 1, \dots, n} (T_{i} Δ_{i})}

, and its respective event indicator is set to be

Δ_{i}^{p} = 0

, as there is no empirical evidence for the potential failure times beyond

\max_{i = 1, \dots, n} {T_{i} Δ_{i}}

When performing a permutation test, the reported p-value can be viewed as an approximation of the true p-value, based on all possible permutations. In the above imputation-based permutation procedure, additional variability in a p-value is expected due to random imputations. To reduce this variability, multiple imputations can be used, such that for each random imputation, B permutations are generated. Assume M imputations are used. Then the p-value is defined as the fraction of the test statistics among the MB test statistics that are at least as large as the observed test statistic Q.

In the following theorem, it is argued that our proposed tests are consistent against all alternatives. The proof is presented in details in the Supplementary Materials.

Theorem:

Let X be a positive failure time random variable, either continuous or discrete, and Y be a categorical random variable with K categories. Let $π_{k} = \lim_{n \to \infty} n_{k} / n, k = 1, \dots, K$ . Assume there are at least two cumulative distribution functions $F_{g} (x)$ and $F_{m} (x), g, m \in {1, \dots, K}$ , such that $F_{g} (x_{0}) \neq F_{m} (x_{0})$ for some $x_{0} \in R^{+}, π_{g} > 0, π_{m} > 0$ , and the conditional censoring distributions are such that $pr (C > x_{0} | Y = g) > 0$ and $pr (C > x_{0} | Y = m) > 0$ . Then, the imputation-based permutation test with the test statistic Q is consistent, namely, the power of the test increases to 1 as $n \to \infty$ .

2.4 Computation time

Table 1 provides the run time of the proposed tests of one dataset, K = 2, under the null hypothesis, one imputation, and 1000 permutations, for different total sample sizes n and $n_{1} = n_{2} = n / 2$ . These results were generated by a 6 years old Intel i7-3770 CPU 3.4 GHz, without paralleling with the different cores of the computer. The first two rows are for identical censoring distributions, and the last two are for different censoring distributions. Evidently, even with n = 1000 observations and low censoring rates, the run time on such a simple computer is not longer than 3.5 min. We expect that an increase of the number of imputations and permutations will increase the run time in a linear fashion.

Table 1.

Computation time (seconds).

Censoring rates of groups 1 and 2	n
Censoring rates of groups 1 and 2	100	200	300	400	1000
(25%, 25%)	1.7	7.1	16.5	30.0	204.2
(50%, 50%)	0.9	3.5	8.0	14.5	97.3
(27%, 55%)	0.9	3.8	8.8	16.1	114.1
(40%, 55%)	0.8	3.2	7.5	13.7	93.8

3 Simulation study

3.1 Simulation design

An extensive numerical study has performed to systematically examine the behavior of our proposed K-sample omnibus non-proportional hazards (KONP) tests under a wide range of two-sided alternatives, various sample sizes, and a wide range of censoring distributions, including unequal censoring distributions. The main part of the simulation study was dedicated to the popular two-sample setting, but settings of K = 3, 4, 5 were considered as well.

As competitors under the 2-sample setting, the following tests were included: the logrank test; Peto–Peto weighted logrank test² that uses a weight function that is very close to the pooled KM estimator; Pepe–Fleming weighted KM test³ with geometric mean of the two KM censoring-distribution estimators as a weight function; and Yang–Prentice test, an adaptive weighted logrank test where the adaptive weights utilize the hazard ratio obtained by fitting the model of Yang and Prentice.⁵ The tests of Uno et al.¹⁶ are invalid under unequal censoring distributions (as demonstrated below), and thus are not included in the following power comparisons.

Table 2 (main text) and Tables 4 and 5 in Appendix 1, provide a comprehensive summary of the 17 non-proportional hazards scenarios and 7 proportional or close to proportional hazards functions, that were studied. For each scenario, the failure and censoring distributions are explicitly provided, and the survival functions of the two groups are plotted. A reference is provided indicating the source of each setting. In short, Scenario A shows differences at mid time points, but similarity in early and late times. Scenarios B–D shows differences in early times. Scenario E is of equal survival functions at early times and of proportional hazards at mid and late times. Scenarios F and G are with crossing hazards. Scenario H is of a U-shape hazards ratio. Scenarios I, J and K are with crossing hazards, based on the following hazards-ratio model of Yang and Prentice⁵

HR (t) = \frac{λ_{2} (t)}{λ_{1} (t)} = \frac{θ_{1} θ_{2}}{θ_{1} + (θ_{2} - θ_{1}) S_{1} (t)}, 0 < t < τ_{0}, θ_{1}, θ_{2} > 0

(1)

where

τ_{0} = \sup {t : S_{1} (t) > 0}

, λ₂ and λ₁ are the hazard functions of the two groups, and S₁ is the survival function of group 1. Under Model (1),

θ_{1} = \lim_{t \to 0} HR (t)

and

θ_{2} = \lim_{t \to τ_{0}} HR (t)

. Also it is assumed that for a continuous function S₁,

HR (t)

is a strong monotone function of t, i.e.

sign {d HR (t) / d t}

is the same for all

t \in (0, τ_{0})

. Scenarios I-1, I-2 and I-3 are with crossing hazards under Model (1). In I-1 the hazard functions cross earlier compared to I-2 and I-3. Scenarios J-1, J-2 and J-3 are with crossing hazards in which the strong monotonicity assumption is violated. In J-1 and J-2 the hazards are piece-wise proportional, and the hazard functions cross in mid time points. In Scenario J-3 the hazard ratio is a continuous function of t, and the hazards cross at a late time point. Under Scenarios K-1, K-2 and K-3 Model (1) is violated, but the strong monotonicity assumption of

HR (t)

holds. In Scenarios K-1 and K-2 the hazards cross at early-mid times, and in Scenario K-3 at mid times.

Table 2.

Unequal censoring distributions.

Group	Small difference	Substantial difference
1	$\min {U (a, b), n">Exp (θ_{1})}$	$\min {U (a, b), n">Exp (θ_{1})}$
2	$\min {U (a, b), n">Exp (θ_{2})}$	U(a, b)

For each scenario described above, four different censoring distributions were considered, two with equal and two with unequal censoring distributions. Under equal censoring distributions, the censoring distributions were taken to be similar to the corresponding referenced paper. Exponential distributions were used for all other scenarios, with ∼25% or 50% censoring rates. Under unequal censoring distributions, the censoring distributions of Wang et al.¹⁵ were used (Table 2). The specific values of $(a, b, θ_{1}, θ_{2})$ are provided in Tables A1 and A2. Under small differences, the censoring rates between the two groups are about 40% and 55%, where 27% and 55% are the respective values under substantial differences. Identical enrollment time was assumed, and the follow-up time ended at event time or censoring time, whichever comes first. The various censoring settings are such that the power of a specific test under a specific scenario is not necessarily increasing as the censoring rate decreases.

Each of the configurations was studied with n = 100, 200, 300 or 400, n₁ = n₂, and performances are summarized based on 2000 replications.

A smaller simulation study was done for K > 2. As competitors, the logrank and Peto-Peto tests were included. For the null scenario K = 3, 4, 5 were studied, and under Scenarios D and J-2, K = 3 was examined. Various sample sizes and a wide range of censoring distributions were considered. A detailed description of these scenarios can be found in the Supplementary Materials.

3.2 The test of Yang and Prentice¹

Since the Yang–Prentice test is the strongest competitor in terms of power, for the to-sample setting, we highlight some of its properties. The Yang–Prentice test is based on Model (1), where the indices 1 and 2 indicate the control and treatment groups, respectively. Since this model is asymmetric in terms of F₁ and F₂, the test is not group-label invariant. By exchanging the group labels between “treatment” and “control,” a different p-value would be provided. This property is unique to this test, and all other tests considered in this work are invariant to group labeling. Consequently, in applications with no clear link between the two groups to treatment/control status (e.g. comparing females versus males, or young versus old), it is unclear how the Yang–Prentice test should be applied.

In order to make the Yang–Prentice test invariant, we apply a permutation test based on the minimum p-value of the two labeling options. Specifically, let PV₁ and PV₂ be the p-values of the Yang–Prentice tests based on $(X_{i}, Δ_{i}, Y_{i}), i = 1, \dots, n$ , and $(X_{i}, Δ_{i}, {\tilde{Y}}_{i}), i = 1, \dots, n$ , respectively, where $Y_{i} \in {1, 2}$ and ${\tilde{Y}}_{i} = I (Y_{i} = 2) + 2 I (Y_{i} = 1)$ . Then, the Yang–Prentice invariant test statistic is defined by $Q_{Y P} = \min {P V_{1}, P V_{2}}$ . The p-value of the imputation-based permutation test, based on the statistic Q_YP, is the fraction of replicates of Q_YP under random permutations of the data, as described in Section 2.3., that are at least as small as the observed test statistic. Since our work mainly focuses on invariant tests, we report and discuss in the main text the performance of the Yang–Prentice invariant test. The results of the original Yang–Prentice test, with the two possible options of labeling for treatment and control, can be found in the Supplementary Materials.

The original Yang–Prentice test (implemented in the R package YPmodel) uses the asymptotic distribution of the test statistic under the null hypothesis. Often, its empirical size is greater than the nominal level, especially under small sample size (see the Supplementary Material). In a recent work of Yang,¹⁷ which deals with interim monitoring using adaptive weighted log-rank test, the author suggests using the re-sampling method of Lin et al.¹⁸ instead of the original asymptotic approach,¹⁸ for improving the type I error rate. The method of Lin et al.¹⁸ is also based on asymptotic results. Alternatively, the imputation approach of Wang et al.¹⁵ works very well in very small sample sizes, and it does not require asymptotic distribution. To make the comparison between the methods consistently, we use the imputation approach for both, our proposed method and the invariant version of Yang–Prentice test.

3.3 Simulation results

Figure 2 provides the empirical power of the tests under the null hypothesis, with equal and unequal censoring distributions. Evidently, under equal censoring distributions all the tests are valid, as the empirical size of the tests are reasonably close to the nominal value 0.05. On the other hand, under the null hypothesis and unequal censoring distributions, the empirical sizes of Uno et al. tests are much higher than the nominal value 0.05. For example, under a sample size of n = 400, and censoring rates of ∼27% and 55%, the empirical size of Uno et al. V₂ bona fide test is 0.099; all the other Uno et al. tests’ sizes are even higher. The empirical sizes of all the other tests are reasonably close to their nominal value. Thus, Uno’s tests will not be considered in the rest of this simulation study.

Figure 2.

Empirical power under the null: top: equal censoring rates; bottom: unequal censoring rates.

Figure 3 summarizes the empirical power of the tests under settings generated by others, while Figure 4 is based on scenarios generated by us. As expected, the power of each test increases with the sample size. In some scenarios the power of the tests increases as the censoring rate increases, since in these settings the non-censored data are centered mainly at the parts of the hazards, which are closer to proportionality and the censored data are mainly located at the non-proportionality area of the hazards. For example, in Scenario I-3 with n = 400, as the censoring rate increases from about 25% to about 50%, the power of Peto–Peto increases by 0.159, Yang–Prentice by 0.2, Pepe–Fleming by 0.218, and logrank by 0.4. The power increase in our KONP tests is much smaller, 0.012 and 0.009.

Figure 3.

Simulation Results of non-proportional hazards settings considered by others: Setting (a) shows differences at mid time points, but similarity in early and late times. (b) to (d) show differences in early times. (e) It is of equal survival functions at early times and of proportional hazards at mid and late times. Scenarios (f) and (g) are with crossing hazards under model (1). Scenario (h) is of a U-shape hazards ratio.

Figure 4.

Simulation results of additional settings of non-proportional hazards: Scenarios I-1, I-2 and I-3 are with crossing hazards under Model (1). In I-1 the hazard functions cross earlier compared to I-2 and I-3. Scenarios J-1, J-2 and J-3 are with crossing hazards in which the strong monotonicity assumption is violated. In J-1 and J-2, the hazards are piece-wise proportional, and the hazard functions cross in mid time points. In Scenario J-3, the hazard ratio is a continuous function of t, and the hazards cross at a late time point. Under Scenarios K-1, K-2 and K-3, Model (1) is violated, but the strong monotonicity assumption of $HR (t)$ holds. In Scenarios K-1 and K-2, the hazards cross at early mid times and in Scenario K-3 at mid times.

Evidently, our KONP tests are often more powerful compared to all other tests. These includes scenarios A, C, D, E, F, J-1, J-2, J-3, I-1, K-1, K-2, K-3. The superiority of KONP over Yang–Prentice under F and I-1 is surprising, since these scenarios follow their Model (1). For example, with n = 400 and 25% censoring rate, the empirical power is about 90% for KONP tests and only 70% for Yang–Prentice. In Scenario G (close to proportional hazards), which also follows Model (1), the results of Yang–Prentice and logrank are similar and often slightly better than KONP tests. In scenarios J-1 and J-2 our tests are substantially more powerful than all the competitors. For example, in Scenario J-2, n = 400, and 25% censoring rate, KONP tests with about 95% power, while the second most powerful test is Yang–Prentice with a power of 48%.

In scenarios G, I-2 and in some of the censoring rates in I-3, the Peto–Peto and Pepe–Fleming tests tend to be with the highest power. In contrast, in many of the scenarios in which the survival functions cross, their power is much lower than our tests. For example, in Scenario K-1, n = 400, and 25% censoring rate, the power of KONP is 89%, while Pepe–Fleming power is 42% and Peto–Peto is 5%.

To conclude, for the 2-sample setting, in most of the non-proportional hazards settings, our proposed KONP tests tend to be more powerful than the other tests, and the differences between S_P and S_LR are very small, if any.

Results of settings with K > 2 can be found in Table S1 and Figure S1 of the Supplementary Material. Based on these results we conclude that all the tests are valid, as the empirical size of the tests are reasonably close to 0.05. For K = 3 and scenarios D and J-2, we see similar results to those shown with K = 2, as KONP tests are often much more powerful than the logrank and Peto–Peto test. For example, for scenario D with n = 200 and 25% censoring rates in all the three groups, the powers of KONP tests are ∼92%, while of Peto–Peto is 49%, and of logrank is 18%.

Figure 5 summarizes the power of the two-sample tests under proportional hazards or close to proportionality. Under these settings, the logrank test is often with the highest power among the invariant tests, as expected. The invariant version of Yang–Prentice test is similar to the logrank test, and the proposed KONP tests, sometimes, have less power.

Figure 5.

Empirical power the two-sample settings under proportional hazards or close to proportionality: Scenario L is with proportional hazards. M is close to proportionality but with substantial differences at early times. N–Q are close to proportionality and model (1).

3.4 A robust approach

Figures 3 and 4 of the main text and Table S2 of the Supplementary Material indicate that under the non-proportional hazards scenarios and among the invariant tests, usually the proposed KONP tests are with the largest power. Under proportional hazards or close to proportionality, usually the logrank test is with the largest power among the invariant tests. The invariant version of Yang–Prentice test is similar to the logrank test under proportional hazards settings, since their model contains the proportional hazards model ( $θ_{1} = θ_{2}$ ).

In case one is interested in a robust powerful test under non-proportional or proportional hazards, the principle of minimum p-value could be adopted based on the elegant Cauchy-combination test of Liu and Xie,¹² which is similar to the test based on the minimum p-value. Denote the p-values of our KONP tests by $p {-value}_{P}$ and $p {-value}_{L R}$ , respectively, and the p-value of the logrank test by $p {-value}_{lgrnk}$ . Then, based on Liu and Xie¹² we define a new test statistic

\begin{array}{l} C a u = [\tan {(0.5 - p {-value}_{P}) π}] / 3 + [\tan {(0.5 - p {-value}_{L R}) π}] / 3 \\ + [\tan {(0.5 - p {-value}_{lgrnk}) π}] / 3 \end{array}

and its p-value is

p {-value}_{Cau} = 0.5 - (\arctan Cau) / π

The candidate tests to be included in Cau are powerful tests under non-proportional (i.e. KONP tests) or proportional hazards (i.e. logrank test and the invariant version of Yang–Prentice test). Due to the similarity in power performances between the logrank test and the invariant version of Yang–Prentice test, under proportional or close to proportional hazards, and due to the high computational burden of Yang–Prentice invariant test, only the logrank test is included.

Figure 6 provides the empirical Type-I error of the two KONP tests, the logrank test and the robust Cau test. Evidently, the size of the tests are reasonably close to 0.05. Figures 7 and 8 summarize the empirical power, based on 1000 replications, of the two KONP tests, the logrank test, the test statistic Cau, along with two additional combined test: the test of Lee¹³ and the MaxCombo test.¹⁴ The test of Lee¹³ is based on the maximum of two weighted logrank test statistics; and the MaxCombo test¹⁴ is based on the maximum of the logrank and three weighted logrank test statistics. (See Supplementary Materials for details.) Table S4 of the Supplementary Materials provides the power values of these tests. Often, the Cau test loses some power comparing to the largest power among KONP, logrank and the other two combined tests, but the loss is relatively small. Lee and MaxCombo tests perform very similarly. In general, our Cau test outperforms Lee and the MaxCombo tests, in terms of power, in all the sub-scenarios (i.e. various sample size and censoring patterns) of type B, C, H, J-1, and J-2; and in some of the sub-scenarios of A, F, I-1, I-2, I-3, K-1, K-2, K-3, and P. All of these are non-proportional hazards scenarios. Under proportional hazards or close to proportionality, i.e. settings L, M, N, O, P, and Q, Lee’s test or MaxCombo outperform Cau. Interestingly, in the cases where Lee or MaxCombo tests are with higher power than Cau, the power loss by using Cau is relatively small; while this is not always the case when Cau outperforms Lee or MaxCombo. For example, under Setting B with n = 200, the power of Cau equals 0.919 while the power values of Lee and MaxCombo are 0.499 and 0.450, respectively.

Figure 6.

Empirical power under the null of the two-sample KONP tests, logrank and the Cau robust test.

Figure 7.

Empirical power under the null of the two-sample KONP tests, logrank, Cau, Lee (2007) and MaxCombo tests: Scenarios A through I-3.

Figure 8.

Empirical power under the null of the two-sample KONP tests, logrank, Cau, Lee (2007) and MaxCombo tests: Scenarios J-1 through Q.

Our R package KONPsurv¹⁹ applies the above robust test Cau as well.

4 Real data examples

4.1 The gastrointestinal tumor data

The Gastrointestinal Tumor Study Group²⁰ compared chemotherapy with combined chemotherapy and radiation therapy, in the treatment of locally unresectable gastric cancer. This dataset was used in Yang and Prentice¹ to demonstrate the utility of their test. Each treatment arm had 45 patients, and two observations of the chemotherapy group and six of the combination group were censored. The primary outcome measure was time to death. The KM survival curves of time to death, in each treatment group are provided in Figure 9. To apply the Yang–Prentice test, we considered chemotherapy as the control group and chemotherapy plus radiation therapy as the treatment group, which is named Yang–Prentice 1. The Yang–Prentice test with reversed group labeling is denoted by Yang–Prentice 2. Table 3 shows the p-values of testing for equality of the survival curves of time to death, of the two treatment groups, against a two-sided alternative, with each of the tests considered in the simulation study. For our tests and the Yang-Prentice invariant test, 10 imputations and 10⁴ permutations for each imputation, were used. Evidently, the smallest p-values are observed under our proposed KONP tests and the Cauchy-combination test Cau.

Figure 9.

Gastrointestinal tumor study: KM curves.

Table 3.

Examples: GST and UCD.

Test	GST p-value	UCD p-value
KONP Pearson	0.0109	0.0049
KONP LR	0.0108	0.0049
Cauchy combination Cau	0.0164	0.0071
Yang–Prentice 1	0.0304	0.0186
Yang–Prentice 2	0.0800	0.0252
Yang–Prentice invariant	0.0479	0.0219
Logrank	0.6350	0.0673
Pepe–Fleming	0.9464	0.2362
Peto–Peto	0.0465	0.3807
Lee	0.0968	0.0054
MaxCombo	0.0908	0.0061

GST: gastrointestinal tumor study; UCD: urothelial carcinoma data.

4.2 Urothelial carcinoma

Few options exist for patients with locally advanced or metastatic urothelial carcinoma after progression with platinum-based chemotherapy. Powles et al.²¹ aimed to assess the safety and efficacy of atezolizumab versus chemotherapy in this patient population. Their study consists of a multi-center, open-label, phase three randomized controlled trial conducted at 217 academic medical centers and community oncology practices mainly in Europe, North America, and the Asia-Pacific region. The primary endpoint was overall survival. Figure S3 of the Supplementary Material of Powles et al.²¹ provides the overall survival KM curves of atezolizumab versus chemotherapy based on 316 and 309 patients, respectively. Although the detailed survival data are unavailable, inspired by Roychoudhury et al.,²² we used Guyot et al.²³ algorithm that maps from digitized curves back to KM data, by finding numerical solutions to the inverted KM equations, using available information on number of events and numbers at risk. The DigitizeIt software was used for reading the coordinates of the KM curves from the published graph. Figure 10 provides the KM curves by treatment arm, and the last column of Table 3 shows the p-values of testing for equality of the survival curves of time to death against a two-sided alternative. Evidently, the smallest p-values are observed under our proposed KONP tests and the second best is the test of Lee. Our Cauchy combined test also performs robustly.

Figure 10.

Urothelial carcinoma data: KM curves.

5 Discussion and conclusions

The proposed KONP tests are based on partition of the sample–space into two subsets, corresponding to three intervals, as of the HHG test.¹⁰ An extensive simulation study shows that when the hazard curves are non-proportional, the KONP tests are often more powerful than all the other tests. In particular, the proposed test is even more powerful than the Yang–Prentice test, under their model with non-proportional hazards. The simulation results show very little differences in power, if any, between the Pearson chi-squared test statistic and the log-likelihood ratio statistic. Since the chi-squared statistic is slightly more powerful, this test statistic is recommended. The proposed Cau test is often more powerful than the other combined tests, such as Lee and MaxCombo tests.

The diverse set of censoring settings considered in this work reflects the complex effect of a censoring distribution on the power: the power of a specific test under a given scenario is not necessarily increasing as the censoring rate decreases.

Other partitions and summary statistics can be easily adopted. In particular, one may consider the extended AD tests of Thas and Ottoy⁹ and Heller et al.¹¹ with higher sample–space partitions and test statistics that aggregate over all partitions by summation or maximization. Nevertheless, for non-censored data, Heller et al.¹¹ showed by simulations (see their Table 1), that increasing the number of partitions can improve power over the HHG 2-sample test under settings in which the density functions intersect 4 times or more. Otherwise, the HHG 2-sample test tends to be more powerful. Figure 11 in Appendix I displays the densities of the 17 non-proportional hazards scenarios studied in this work. Evidently, the survival scenarios considered by others and by us are of less than four intersections. Simulation study of the AD test statistic with sample–space partition of two intervals, yields lower power than the proposed KONP tests. A comprehensive comparison with other sample–space partitions and aggregations could be a topic of future research.

This work suggests tests that accommodate right-censored data. Since the tests are based on the KM estimator, it seems that a modification to left truncation might be possible. However, additional work is required to modify the imputation-permutation approach for left-truncated data.

Implementation of our tests, KONP-P, KONP-LR, and Cau, is available in the R package KONPsurv,¹⁹ which can be freely downloaded from CRAN.

Supplemental Material

SMM907355 Supplemental Material - Supplemental material for K-sample omnibus non-proportional hazards tests based on right-censored data

Supplemental material, SMM907355 Supplemental Material for K-sample omnibus non-proportional hazards tests based on right-censored data by Malka Gorfine, Matan Schlesinger and Li Hsu in Statistical Methods in Medical Research

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: NIH (R01CA189532) and the U.S.-Israel Binational Science Foundation (2016126).

ORCID iD

Malka Gorfine

Supplementary material

Supplementary material available at SMMR online includes: (1) The proof of the theorem. (2) Description of some of the tests included in the simulation study. (3) Description and results of K-sample simulation settings with K = 3, 4, 5. (4) Plots of the empirical power results of the robust test Cau. (5) The exact empirical power used for generating all the plots in the main text and the Supplementary Material file.

Appendix 1. Detailed description of the simulation settings

References

Yang

Prentice

RL.

Improved logrank-type tests for survival data using adaptive weights.

Biometrics 2010; 66: 30–38.

Peto

Asymptotically efficient rank invariant test procedures. J R Stat Soc A 1972; 185–207.

Pepe

Fleming

TR.

Weighted Kaplan–Meier statistics: a class of distance tests for censored survival data. Biometrics 1989; 497–507.

Pepe

Fleming

TR.

Weighted Kaplan–Meier statistics: large sample and optimality considerations. J R Stat Soc B 1991; 341–352.

Yang

Prentice

RL.

Semiparametric analysis of short-term and long-term hazard ratios with two-sample survival data. Biometrika 2005; 92: 1–17.

Darling

DA.

The Kolmogorov–Smirnov, Cramer–von Mises tests. Ann Math Stat 1957; 28: 823–838.

Pettitt

AN.

A two-sample Anderson–Darling rank statistic. Biometrika 1976; 63: 161–168.

Scholz

Stephens

MA.

k-sample Anderson–Darling tests. J Am Stat Assoc 1987; 82: 918–924.

Thas

Ottoy

JP.

An extension of the Anderson–Darling k-sample test to arbitrary sample space partition sizes. J Stat Comput Simul 2004; 74: 651–665.

10.

Heller

Gorfine

A consistent multivariate test of association based on ranks of distances. Biometrika 2013; 100: 503–510.

11.

Heller

Kaufman

, et al. Consistent distribution-free K-sample and independence tests for univariate random variables. J Mach Learn Res 2016; 17: 978–1031.

12.

Liu

Xie

Cauchy combination test: a powerful test with analytic p-value calculation under arbitrary dependency structures. J Am Stat Assoc 2019; 1–18.

13.

Seung-Hwan

On the versatility of the combination of the weighted log-rank statistics. Comput Stat Data Anal 2007; 51: 6557–6564.

14.

Lin

Roychoudhury

, et al. Alternative analysis methods for time to event endpoints under non-proportional hazards: a comparative analysis. Statistics in Biopharmaceutical Research 2020; pp.1–19.

15.

Wang

Lagakos

Gray

RJ.

Testing and interval estimation for two-sample survival comparisons with small sample sizes and unequal censoring.

Biostatistics 2010; 11: 676–692.

16.

Uno

Tian

Claggett

, et al.

A versatile test for equality of two survival functions based on weighted differences of Kaplan–Meier curves.

Stat Med 2015; 34: 3680–3695.

17.

Yang

Interim monitoring using the adaptively weighted log-rank test in clinical trials for survival outcomes.

Stat Med 2019; 38: 601–612.

18.

Lin

Wei

Ying

Checking the cox model with cumulative sums of martingale-based residuals. Biometrika 1993; 80: 557–572.

19.

Schlesinger

Gorfine

Konp tests: powerful k-sample tests for right-censored data. R package version 1.0.1. 2019.

20.

Schein

and Gastrointestinal Tumor Study Group. A comparison of combination chemotherapy and combined modality therapy for locally advanced gastric carcinoma. Cancer 1982; 49: 1771–1777.

21.

Powles

Durán

van der Heijden

, et al.

Atezolizumab versus chemotherapy in patients with platinum-treated locally advanced or metastatic urothelial carcinoma (IMvigor211): a multicentre, open-label, phase 3 randomised controlled trial.

Lancet 2018; 391: 748–757.

22.

Roychoudhury

Anderson

, et al. Robust design and analysis of clinical trials with non-proportional hazards: a straw man guidance from a cross-pharma working group. arXiv preprint 2019; DOI: arXiv:1908.07112.

23.

Guyot

Ades

Ouwens

, et al.

Enhanced secondary analysis of survival data: reconstructing the data from published Kaplan-Meier survival curves.

BMC Med Res Methodol 2012; 12(1): 9.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.43 MB

K -sample omnibus non-proportional hazards tests based on right-censored data

Abstract

Keywords

1 Introduction

2 K-sample tests based on sample–space partition

2.1 Motivation and notation

2.2 The test statistic

2.3 The permutation procedure

2.4 Computation time

3 Simulation study

3.1 Simulation design

3.2 The test of Yang and Prentice 1

3.3 Simulation results

3.4 A robust approach

4 Real data examples

4.1 The gastrointestinal tumor data

4.2 Urothelial carcinoma

5 Discussion and conclusions

Supplemental Material

SMM907355 Supplemental Material - Supplemental material for K-sample omnibus non-proportional hazards tests based on right-censored data

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

Supplementary material

Appendix 1. Detailed description of the simulation settings

References

Supplementary Material

3.2 The test of Yang and Prentice¹