Comparison of methods to testing for differential treatment effect under non-proportional hazards data

1.   Introduction

The most commonly used two-sample nonparametric test statistic is the log-rank (LR) statistic, which is locally most powerful against proportional hazards (PH) alternatives. However, rank-based log-rank statistics may not be sensitive to certain types of stochastic ordering alternatives, particularly if the hazard functions cross or a delayed effect is observed. In reality, the PH does often not hold true as for example in some immuno-oncology therapies ^[1]. As an alternative to generalized linear rank statistics, ^[2] proposed nonparametric weighted Kaplan-Meier (WKM) statistics, which are non-rank-based statistics and therefore sensitive to the magnitude and duration of the difference in survival curves over time. These test statistics seems to compare favorably with the popular log-rank test statistic across a broad range of stochastic ordering alternatives.

As a particular case of WKM, you get the restricted mean survival time (RMST) which is an alternative robust and clinically interpretable summary measure that does not rely on the proportional hazard (PH) assumption. ^[3] concluded that the RMST-based test has better performance than the log-rank test when the truncation time is reasonably close to the tail of the observed curves under non-PH scenarios where late separation of survival curves is observed.

Under non-proportional alternatives, the effect of a treatment may wane over time, leading to a decreasing hazard ratio and a closing up of the two survival curves or may have a delayed effect whereby they do not separate until a certain interval of time has elapsed. As another more powerful alternative over the unweighted log-rank statistic, ^[4] proposed a class of adaptive weighted log-rank statistics. These tests may assign more weight on either early, middle or late survival differences with a suitable choice of weight function. A particularly useful family in the class of weighted log-rank statistics was introduced in ^[5] and ^[6]. However, in many situations, it is difficult to know in advance the kind of survival differences that will actually occur. Therefore, ^[7] proposed the maximum of a finite collection of these weighted log-rank statistics trying to overcome the above problem.

Later, ^[8,9] proposed tests based on a more extended family of test statistics where the weights are governed by two parameters and combinations of them that are more sensitive to nonproportional hazard alternatives. Furthermore, ^[10] considered a new combination of these test statistics. On the other hand, ^[11] considered the same family of weighted functions for the WKM test statistic and evaluated the maximum of nine of these test statistics. ^[12] proposed a different weighted functions for these test statistics. Recently, ^[13] investigated weighted log-rank tests, the supremum log-rank test and composite tests (mainly based on log-rank tests). On the other hand, ^[14] studied the performance of nine tests. Seven out of nine are based on weighted log-rank and the other two are WKM and RMST. Both papers came to the conclusion that: there is not a single most powerful test across all scenarios but the maximum of several weighted log-rank tests seems to be the best choice. The former proposed the so-called MaxCombo test and the later the proposal studied by ^[10]. Recently, ^[15] have performed additional simulations to evaluate the MaxCombo test under some extreme scenarios which are unlikely in real life. Futhermore, they have provided design and analysis considerations based on a combination test under different non-PH types and present a straw man proposal for practitioners.

To our knowledge, the class of maximum weighted log-rank tests and that of maximum weighted Kaplan-Meier tests have not been compared to each other. Therefore, we try to fill this gap as well as to compare them with the RMST which is easily interpretable. The paper is organized as follows: In Section 2, we describe the tests to be compared and our approach to get the critical points. A simulation study of the performance of the tests is presented in Section 3. In Section 4, we apply all the methods to a real data set from a clinical trial of marrow transplant patients with leukemia. We finish with a discussion in Section 5.

2.   Methods

We assume the two-sample general random censorship model. Consider two censored samples with sample sizes $n_{1}$ from the treatment group and $n_{2}$ from the control group. Let $T_{ij},$ $i = 1, 2,$ $j = 1, ..., n_{i},$ be independent, positive random variables and $C_{ij}$ be independent censoring variables that are also independent of the survival variables $T_{ij}.$ The data really observed are $\left(X_{ij}, \delta _{ij}\right),$ $i = 1, 2,$ $j = 1, ..., n_{i},$ where $X_{ij} = \min (T_{ij}, C_{ij})$ and $\delta _{ij} = I\left(T_{ij}\leq C_{ij}\right),$ where $I$ is the indicator function. Define survival functions $S_{i}(t) = P(T_{ij}\geq t)$ and $C_{i}(t) = P(C_{ij}\geq t), i = 1, 2,$ which are assumed to be absolutely continuous throughout the paper. The null hypothesis to test is that the survival distributions between two groups are the same, that is,

2.1.   Weigthed Log-rank tests and combinations

The so-called weighted log-rank statistics are given by

where $\overline{N}_{i}(t)$ is the number of failures in group $i$ before or at time $t$ and $\overline{Y}_{i}(t)$ the number at risk in group $i$ at time $t-,$ $i = 1, 2, \widehat{w}(t)$ is a bounded nonnegative weight function and $t_c = \sup\{t/\hat C_1(t) \wedge \hat C_2(t) > 0\}$ where $\widehat{C}_{i}$ denotes the Kaplan-Meier estimator of the censoring survival function in group $i$. As the weight function is sensitive to the alternative hypothesis, it should be chosen carefully.

^[7] proposed

with $\widehat{S}(t-)$ being the left-continuous version of the product-limit estimator (^[16]) for the survival function based on the pooled survival data. We consider the family of test statistic given in (2.2) when the class of functions (2.3) are chosen, WLR$^{\lambda, \gamma}$. In particular, $\lambda = \gamma = 0$ corresponds to the log-rank statistic, $WLR^{0, 0}$, which specifies equal weights and $\lambda = 1$ and $\gamma = 0$ corresponds to the Prentice-Wilcoxon statistic, $WLR^{1, 0}$, which places more weight on the earlier time points. In contrast, $WLR^{0, 1}$ places more weight on the later time points and $WLR^{1, 1}$ will emphasize the middle differences. ^[8] studied the maximum over the statistics $WLR^{0, 0}, WLR^{2, 0},$ $WLR^{0, 2}$ and $WLR^{2, 2}$ as well as their average. ^[9] studied the test statistics $\left\vert WLR^{1, 0}+WLR^{0, 1}\right\vert /2,$ $\left(\left\vert WLR^{1, 0}\right\vert +\left\vert WLR^{0, 1}\right\vert \right) /2$ and $\max \left(\left\vert WLR^{1, 0}\right\vert, \left\vert WLR^{0, 1}\right\vert \right)$. ^[10] and ^[13] considered $WLR_{max3} = \max \left(\left\vert WLR^{0, 0}\right\vert, \left\vert WLR^{1, 0}\right\vert, \left\vert WLR^{0, 1}\right\vert \right)$. ^[14] focused on the so called MaxCombo test $WLR_{max4} = \max (\left\vert WLR^{0, 0}\right\vert,$ $\left\vert WLR^{0, 1}\right\vert, \left\vert WLR^{1, 0}\right\vert, \left\vert WLR^{1, 1}\right\vert)$.

According to the conclusions of the papers that have studied these test statistics and combinations, we focus on maximum combination tests. To conduct inference on these test statistics, i.e. to calculate the p-value for testing (2.1), there are different methods. One is to obtain the asymptotic distribution of these test statistics as the maximum of a multivariate normal distribution with mean zero and covariance that can be consistently estimated by

for $i, j = 1, 2, 3, 4, \hat w_i(t), i = 1, 2, 3, 4$ are the weight functions with $\lambda = \gamma = 0; \lambda = 1, \gamma = 0; \lambda = 0, \gamma = 1$ and $\lambda = \gamma = 1$, respectively and $\Delta \overline{N}_{i}(t) = \overline{N}_{i}(t)- \overline{N}_{i}(t-)$, the number of events at time $t$ in group $i$. However, in this paper, to account for small and moderate sample sizes, bootstrap is used to obtain p-values as is explained in Section 2.3.

2.2.   Weigthed Kaplan-Meier tests and combinations

^[2] defined the weighted Kaplan-Meier statistics as

They proposed using $\widehat{u}(t) = \widehat{C}_{1}(t-)\widehat{C} _{2}(t-)/(p_{1}\widehat{C}_{1}(t-)+p_{2}\widehat{C}_{2}(t-))$ to stabilize the test statistic and down weight its variance toward the end of the observation period if censoring is heavy.

^[11] proposed

with $p_{i} = n_{i}/(n_1+n_2)$ and $\widehat{C}_{i}$ denotes the Kaplan-Meier estimator of the censoring random variables in group $i$ and $\hat w(t)$ is defined in (2.3). We shall call these test statistics $WKM^{\lambda, \gamma}$. In particular, $\lambda = \gamma = 0,$ $WKM^{0, 0},$ corresponds to the Pepe-Fleming's WKM test. ^[11] studied $\max_{i, j = 0, 1, 2}\left(WKM^{i, j} \right).$

As a particular case, we have the difference in RMST (^[3])

To be fair the comparison of these test statistics with that based on the log-rank test, we will focus on their maximum combination tests counterparts. In particular, $WKM_{max4} = \max (WKM^{0, 0}, WKM^{0, 1},$ $WKM^{1, 0}, WKM^{1, 1})$ and $WKM_{max3} = \max (WKM^{0, 0}, WKM^{0, 1}, WKM^{1, 0})$. To conduct inference on these test statistics, i.e. to calculate the p-value for testing (2.1), the asymptotic distribution of these test statistics can be obtained as the maximum of a multivariate normal distribution with mean zero and covariance that can be consistently estimated by

for $i, j = 1, 2, 3, 4$, and $\hat u_{i}(t), i = 1, 2, 3, 4$ are the weight functions given in (2.4) for $\hat w_i(t)$ defined in the previous section. However, in this paper, to account for small and moderate sample sizes, bootstrap is used to obtain p-values as is explained in Section 2.3.

2.3.   Computing p-values using bootstrap

The p-value of our test statistics can be obtained using their asymptotic distributions either by numerical integration or Monte Carlo estimation of the multivariate Gaussian distribution. ^[11] pointed out the variance-covariance matrix may not have a close form solution and requires intense computation. Therefore, our proposal is to use Bootstrap instead of asymptotic distributions that has less computational cost and improve the precision of the asymptotic approximations in small and moderate samples as well as deal with analytically challenging problems in some cases.

We propose to calculate the p-values as follows:

● Step 1: To obtain the value of the test statistic for the original sample, $T_0$

● Step 2: Draw B bootstrap samples from the original data, so $n_1+n_2$ observations are obtained and consider the first $n_1$ as controls and the other $n_2$ as treatments then to obtain the value of the test statistic for the bootstrap sample, $T_B(i), i = 1, ..., B$.

● Step 3: Estimate the p-value by counting the number of $T_B(i)$ greater than $T_0$.

We also offer R-code facilitating the implementation of the evaluation of the p-value of the test statistics in Supplementary.

3.   Simulation study

We carry out a simulation study to compare the performance of these thirteen test statistics: $WLR^{0, 0}$, $WLR^{0, 1}$, $WLR^{1, 0}$, $WLR^{1, 1}$, $WLR_{max4}$, $WLR_{max3}$, $WKM^{0, 0}$, $WKM^{0, 1}$, $WKM^{1, 0}$, $WKM^{1, 1}$, $WKM_{max4}$, $WKM_{max3}$ and $D$. The simulation design described in ^[9] has been considered but it is similar to that included in the other cited references such as ^[17] and ^[13]. We used piecewise exponential models to generate simulated data with parameters $\lambda_i(t)$ from groups 1 and 2 set to represent common alternatives that might arise in real data. To be specific, in Figure 1, (a) null case $\lambda_1 = 2, \lambda_2 = 2$, (b) proportional hazards $\lambda_1 = 1, \lambda_2 = 2$, (c) early survival differences, $t < 0.3$, $\lambda_1 = 3, \lambda_2 = 0.75$; $0.3 \leq t < 0.6$, $\lambda_1 = 0.75, \lambda_2 = 3$; $t \geq 0.6$, $\lambda_1 = 1, \lambda_2 = 1$, (d) late survival differences, $t < 0.5$, $\lambda_1 = 2, \lambda_2 = 2$; $t \geq 0.5$, $\lambda_1 = 4, \lambda_2 = 0.4$ and (e) early and late occurring survival differences, $t < 0.2$, $\lambda_1 = 3, \lambda_2 = 0.75$; $0.2 \leq t < 0.4$, $\lambda_1 = 0.75, \lambda_2 = 3$; $0.4 \leq t < 0.6$, $\lambda_1 = 1, \lambda_2 = 1$; $t \geq 0.6$, $\lambda_1 = 3, \lambda_2 = 0.75$.

In each case, the censoring distributions are uniform and the censoring percentages in the sample are $0%, 20%, 40%$ and $60%$. Empirical size and power of the tests were evaluated considering $\alpha = 0.05$ for sample sizes $n_1 = n_2 = 20, 50, 70,100$ and $150$. For each sample size, 2000 replications were performed for each configuration of survival and censoring distributions and in each one of them 500 bootstrap repetitions were performed.

The simulations were conducted using R (https://www.r-project.org//R Core Team, Vienna, Austria), specifically, we simulate the data with the $\mathtt{rpwexp}$ function of the $\mathtt{nhpsim}$ package ^[18], which simulates the exponential distribution by parts and allows to specify a distribution where the risk rate changes with time, $\mathtt{logrank}$.$\mathtt{test}$ and $\mathtt{logrank}$.$\mathtt{maxtest}$ functions of the $\mathtt{nph}$ package ^[19] were used to calculate the statistics $WLR^{0, 0}$, $WLR^{0, 1}$, $WLR^{1, 0}$, $WLR^{1, 1}$, $WLR_{max4}$, $WLR_{max3}$ and the $\mathtt{survtest}$ function of the $\mathtt{SurvBin}$ package ^[20] to calculate the statistics $WKM^{0, 0}$, $WKM^{0, 1}$, $WKM^{1, 0}$, $WKM^{1, 1}$ and we implement a function to calculate the maximum weighted Kaplan-Meier test statistics, $WKM_{max4}$ and $WKM_{max3}$. Finally, we calculate the D statistic with the $\mathtt{rmst2}$ function of the $\mathtt{survRM2}$ package ^[21].

As can be seen in Table 1, the simulated type Ⅰ error of the thirteen test statistics for all sample sizes are close to the nominal, two-sided 5% significance level. Therefore, we can assess their power performance.

The power of the thirteen test statistics for each simulation scenario is shown in Table 2. Figure 2 displays the powers under 60% censoring. Under proportional hazard configuration, the log-rank test ($WLR^{0, 0}$) and Pepe-Flemings's weighted Kaplan-Meier test ($WKM^{0, 0}$) have maximum power as expected. However, the power of $WKM_{max3}$ and $WKM_{max4}$ are really close. The D statistic is also competitive but under non-censoring. Therefore, these two combination tests have almost the same sensitivity as $WLR^{0, 0}$, i.e. the optimal test against proportional hazards alternative.

For early survival differences, as expected, $WLR^{1, 0}$ and $WKM^{1, 0}$ have greater power than the other statistics. Even though, $WLR_{max3}$ fares poorer than them, it is the best among the rest of considered test statistics. In the case of the late survival differences, $WLR^{0, 1}$ and $WKM^{0, 1}$ have greater power than the others. Nevertheless, the difference in power between these and $WKM_{max3}$ is small. Finally, in the case of early and late survival differences, there is not a clear more powerful test statistic but $WKM_{max3}$ is consistently a good choice for all censoring and sample sizes.

In summary, assuming ignorance of the likely treatment effect, except for early survival differences, the $WKM_{max3}$ test statistic is preferable to the other maximum combination test statistics. For early survival differences, the $WLR_{max3}$ test statistic exhibits a higher power than $WKM_{max3}$ but the latter is higher than the well-known log-rank test. In fact, ^[14] showed a minimal power gain power for WKM and RSMT to the LR test when there is a delayed effect but $WKM_{max3}$ has better performance than WKM and RMST.

4.   Application to a real data set

We ilustrate the performance of the thirteen test statistics with a real data set, chosen because they appear to show late difference. The data is available in the package $\mathtt{KMsurv}$ ^[22] with the name $\mathtt{bmt}$. The $\mathtt{bmt}$ data come from a study on the survival of bone marrow transplant patients with different types leukemia. The complete data frame has 22 variables and to carry out the study we worked with the following: disease free survival time, that is, the time either relapse or death (t2), censoring indicator, 1: relapse or death, 0: censored survival time (d3), disease group, acute lymphoblastic leukemia (ALL) consisting of 38 patients and low risk acute myeloid leukemia (AML) consisting of 54 patients.

Kaplan-Meier curves are shown in Figure 3 and exhibit a late difference in survival. The p-values associated to the proposed test statistics are summarized in Table 3, from which it can be seen that most of the test statistics reject the null hypothesis $H_0: S_1 = S_2$ at the 0.05% level. However, $D$, $WLR^{0, 1}$, $WLR_{max3}$ and $WLR_{max4}$ yield a nonsignificant p-value. Although, $WLR_{max3}$ and $WLR_{max4}$ p-values are in the borderline. Therefore, $WKM_{max3}$ and $WKM_{max4}$ yield stronger evidence against the null hypothesis of no difference between the two groups than $WLR_{max3}$ and $WLR_{max4}$. Therefore, one can miss a potential regulatory submission opportunity if the primary analysis is carried out based on the well-known RMST or MaxCombo tests ($WLR_{max4}$ or $WLR_{max3}$).

5.   Discussion

In practice, the log-rank test is the most used to test the equality of two survival distributions in the presence of censoring and it is known to have optimum power to detect a difference in survival distributions when proportional hazards holds. However, non-proportional hazard ratio between the control and experimental arms occurs fairly often in clinical trials. For example, the proportional hazards (PH) assumption often does not hold for the primary time-to-event (TTE) efficacy endpoint in trials of novel immuno-oncology drugs. For some immuno-oncology drugs, a delayed treatment effect is observed, where the survival curves for the immuno-oncology drug and standard of care are not separated initially but start separating after some period of time. In other oncology trials, the separation between the survival curves occurs early on but then the distance between the curves diminishes over time. There may also be cases where the survival curves cross each other, i.e., the PH assumption is violated. With a preconceived type of treatment effect, then the weighted log rank and weighted Kaplan-Meier test statistics are good alternatives to the log-rank test. However, if the direction of the violation from non-PH is not known, then a combination of the previous test statistics should be used. There are a number of papers studying these combinations based on weighted log rank and weighted Kaplan-Meier test statistics separately. Recently, ^[13] and ^[14] used Monte Carlo simulation to study the performance of nine different statistics based mainly on weighted log-rank tests. They conclude that a modified version of $WLR_{max3}$ and $WLR_{max4}$, respectively are preferred in the absence of prior knowledge regarding the PH or non-PH patterns. ^[15] concluded that MaxCombo test shows clear advantages when a large treatment benefit emerges later on the trial. They introduced a design approach for confirmatory trials with $WLR_{max4}$ and developed a stepwise and iterative approach for calculating sample size when the final analysis is based on this test statistic. Our conclusion is that $WLR_{max3}$ is better than the MaxCombo test, $WLR_{max4}$ and $WKM_{max3}$ is preferable to $WLR_{max3}$ in most of the cases for small and moderate sample sizes. Performance of the KM based methods is known that depend on the lenght of study period ($t_c$) and censoring pattern. We have minimum of the largest observed time in each of the two groups as a choice of $t_c$. There are other issues such as the interpretation of the WLR test statistics result in terms of the clinical effect. ^[23] proposed a Cox-model based time-varying treatment effect estimate to complement the WLR test statistics. However, investigation of these topics is beyond the scope of the present paper.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The second author was supported by grant PID2019-106861RB-I00 funded by MCIN/AEI/10.13039/501100011033 and grant A-SEJ-154-UGR20 funded by FEDER/Junta de Andalucía/Department of Economic Transformation, Industry, Knowledge and Universities and the first author was supported by grant PID2019-104681RB-I00 funded by MCIN/AEI/10.13039/501100011033.

Conflict of interest

The authors declare there is no conflict of interest.

Supplementary

R-code for the p-value calculation

$\mathtt{library(KMsurv)}$

$\mathtt{data(bmt)}$

$\mathtt{B=500}$

$\mathtt{bmt2=datos[(bmt$ group==1)|(bmt $group==2),]}$

$\mathtt{df=data.frame(time=bmt2$ t2, cens=1-bmt2 $d3,}$

$\mathtt{event=ifelse(bmt2$ d3==1, TRUE, FALSE), group=bmt2 $group)}$

$\mathtt{n=nrow(df)}$

$\mathtt{ng1=as.numeric(table(df$ group)[1])}$

$\mathtt{ng2=as.numeric(table(df $group)[2])}$

$\mathtt{max\_wlr=logrank\_maxtest(df$ time, df $event, df$ group,}$

$\mathtt{rho = c(0,0,1,1), gamma = c(0,1,0,1))}$

$\mathtt{Test\_wlr00=max\_wlr $tests$ z[1]}$

$\mathtt{est\_wlr00=vector()}$

$\mathtt{for (b \;in \;1:B)\{}$

  $\mathtt{index=sample(1:n, replace=FALSE)}$

  $\mathtt{dfindex=df[index,]}$

  $\mathtt{g1=dfindex[1:ng1,]}$

  $\mathtt{g1 $group=1}$

  $\mathtt{g2=dfindex[(ng1+1):(ng1+ng2),]}$

  $\mathtt{g2$ group=2}$

  $\mathtt{combinedB=rbind(g1,g2)}$

  $\mathtt{max\_wlr=logrank\_maxtest(combinedB $time, combinedB$ event,}$

  $\mathtt{combinedB $group, rho = c(0,0,1,1), gamma = c(0,1,0,1))}$

  $\mathtt{est\_wlr00[b]=max\_wlr$ tests $z[1] }$

$\mathtt{\}}$

$\mathtt{for(j \;in \;1:B)\{}$

  $\mathtt{if(abs(est\_wlr00[j])>abs(Test\_wlr00)) n\_wlr=n\_wlr+1}$

$\mathtt{\}}$

$\mathtt{prop\_wlr=n\_wlr/ B }$