From September 1, 1991 to December 31, 2005, 758 lung cancer patients underwent tumor removal surgery at Samsung Medical Center in Korea. This study was conducted over more than 16 years and the observed follow-up times ranged from 0.3 to 195 months. In this study, relapse after the surgical process corresponds to the event of interest, and death corresponds to a CR. Among the patients, 580 relapsed (76.5%), 65 died without relapse (8.6%), and 113 were censored (14.9%). The prognostic factor was tumor size (unit: cm) at surgery, which ranged from 0 to 19. A box plot of tumor size with the five numbers is displayed in Fig. 1. Moreover, summary statistics of time to event occurrence and tumor size are presented in Table 1 depending on the event type, such as censoring, relapse, or death.

Let X, Y, and C be the time of the event of interest, the time of CR, and a censoring time, respectively. Let δ=1 if the event of interest occurred, 2 if the CR occurred, and 0 if censored. Let

be observed data, where Z_i is a continuous prognostic factor. Suppose that

are the ordered distinct times at which the event of interest occurs and that

is the set of labels of subjects that fail at t_(u), u=1, 2, …, m. In addition, suppose that

is the set of labels of subjects censored or failed due to CR in the interval [t_(u), t_(u+1)), u=0, 1, …, m.

For a fixed value of µ, set g=1 if Z≤µ and 1 otherwise. Let d_g,u(µ) and r_g,u(µ) denote the number of subjects to whom the event of interest occurs at t_(u) and the number of subjects who are at risk up to time t_(u) in group g, respectively. Let r_u=r_1,u(µ)+r_2,u(µ). Let F_g,µ(t) be the cumulative incidence function (CIF) for the event of interest and S_g,µ(t) be the survival function of being free from any event at t in group g, respectively. The estimated CIF⁸ of F_g,µ(t) is given by

where Ŝ_g,µ(t) is the Kaplan-Meier estimator for S_g,µ(t).10 Define

as a correction factor at t_(u) in group g. Let r̃_g,u(µ)=w_g,u(µ)r_g,u(µ) and r̃_u(µ)=r̃_1,u(µ)+r̃_2,u(µ). Define the score for a subject in D(t_(u)) at t_(u) as

where I_i,g(µ)=1 if subject i is a member of group g and 0 otherwise. Also, define the score for a subject in C(t_(u)) at t_(u) as

for l′∈C(t_(u)). We call these scores the Gray-statistic-type scores and denote ā_µ by their average score.

The testing problem of interest is the independence of X and Z, setting up the null hypothesis as

for all t∈(0,∞) and all cutpoints µ in the prognostic factor Z. It was showed that, for a fixed value of µ, the Gray's test statistic78 for testing the hypothesis of (1) is equivalent to a linear rank statistic,6

Let $H_z\left(\mu \right)=n^{-1}{\sum }_{i=1}^{n}I\left\{Z_i\le \mu \right\}$ denote the empirical distribution function of the prognostic factor Z. Let m_µ=nH_Z(µ) and n_µ=n−m_µ. Following the arguments,125 under H₀, given the scores, the expectation and the variance of L_µ are

and

Then, the standardized statistic of L_µ is defined as

Under H₀, T_µ is asymptotically normally distributed using the arguments.11

Let ξ(ε,Z)=min {µ:H_z(µ)≥ε}. We restrict the possible cutpoints to an interval [ξ(ε₁,Z),ξ(ε₂,Z)] with 0<ε₁<ε₂<1. In the same fashion,125 we define a cutpoint estimator µ̂(ε₁,ε₂) as the value of µ that yields the maximum of the absolute of T_µ, i.e.,

In addition, we define a maximally selected rank statistic as

From the arguments,125 the asymptotic distribution of Q(ε₁,ε₂) is equivalent to the distribution of the supremum of a standardized Brownian bridge. Thus, the approximation12 under H₀ gives

where ϕ(·) denotes a standard normal density.

However, as shown in Table 2, the testing procedure based on the approximation of (2) seems to be conservative. Instead, we propose a permutation test based on B permuted samples by permutingthe observed values of the prognostic factor Z. To be specific, denote Q_b(ε₁,ε₂) by the b^th(b=1, 2, …, B) permuted value of Q(ε₁,ε₂) based on

where P_n is a set of all permutations of integers 1 to n. Suppose that the value of Q(ε₁,ε₂) based on the observed data set {D_i:i=1, 2, …, n} is q₀. Then, the probability in the left-hand side of (2) is empirically determined by

called the empirical p value corresponding to the observed value of q₀.

We performed simulations to investigate the finite-sample performance of the proposed methods in terms of bias and standard deviation (SD) of the cutpoint estimator and the empirical power of the maximally selected test statistics. We generated a time to the event of interest (X) and time to a CR (Y) using Gumbel's bivariate exponential distribution1314 with degree of dependency α set as 0 and 0.3. We generated a non-informative censoring time C from an exponential distribution with a hazard rate of λ(>0), which was determined to satisfy P(C<X⋀Y)=p, where p denotes the censoring fraction. We set p=0 and 0.3. We also generated the prognostic factor Z from a uniform distribution U(0,1) and set the true cutpoint value µ as 0.5. We set ε₁=0.1 and ε₂=0.9. The effect size θ=exp(β) was the relative risk between the two groups of patients, where Z was greater than or less than or equal to µ. We set β as 0, 1, 1.5, 2, and 3: β=0 corresponded to the null hypothesis and β=1, 1.5, 2, and 3 corresponded to the alternative hypotheses. We considered four combinations of α and p: (α,p)=(0,0), (0,0.3), (0.3,0), (0.3,0.3). We performed 500 replications for each configuration of α and p with sample sizes of 50, 100, and 200. We also permuted each sample 250 times to obtain the empirical null distribution of Q(ε₁,ε₂). Fig. 2 depicts the empirical p values of the proposed test Q(ε₁,ε₂) based on a simulated sample against the number b of permutation times for each combination of α and p under H₀:β=0. As shown in Fig. 2, the number B=250 of permutation times was chosen as acceptable regardless of the sample size and the combinations of α and p. See the work6 for details of the data generation procedures. Fig. 3 depicts the empirical distribution function of μ̂ based on 2000 replications, with the sample size of 100 for each combination of α and p when β=0, 1, 1.5, 2, and 3. As expected, regardless of the combinations of α and p, the distribution of μ̂ is centered around the true cutpoint value 0.5 of µ as β increases from 0 to 3.

Table 2 displays the proportion of CR, the bias (Bias) and SD of the cutpoint estimator μ̂, and the approximated (A.Pow) and permutation-based (P.Pow) power of the proposed test statistic Q(ε₁,ε₂). These are based on 500 replications and 250 permutations for each combination of α and p when n=50, 100, and 200. As expected, as β increases, the SD of μ̂ decreases gradually, and both A.Pow and P.Pow converge to 1 regardless of the combination of α and p and sample size n. For a fixed value of β, as n increases, both Bias and the SD of μ̂ decrease; also, both A.Pow and P.Pow increase for any combination of α and p. The permutation-based test (P.Pow), included in the interval (0.031,0.069), satisfies a significance level of 0.05 for most cases, while the approximated test (A.Pow) is very conservative. Furthermore, P.Pow is larger than A.Pow regardless of the combination of α and p and sample size n.

Based on the proposed test statistic Q(ε₁,ε₂), we split the patients into two groups for relapse, a high-risk group and a low-risk group, to apply different treatments to each group. We set ε₁ and ε₂ as 0.1 and 0.9, respectively. As shown in Table 3, at the first split of 758 patients, the criterion was tumor size of 1.8, and the CIFs for relapse of the two groups were different (p value<0.001). We further investigated whether the patients within each subgroup were homogeneous in experiencing relapse after surgery in terms of Q(ε₁,ε₂). The group with a tumor size less than or equal to 1.8 was statistically homogeneous (p value=0.096), while the group having tumor size greater than 1.8 was not homogeneous (p value<0.001). In the same way, we applied a binary split of only the group with a tumor size greater than 1.8. The next split criterion was a tumor size of 3, and only the group with a tumor size greater than 3 was not homogeneous (p value=0.020). The third split was 6.5, and both groups, tumor size less than or equal to 6.5 or greater than 6.5, were homogeneous with p values of 0.530 and 0.957, respectively. The top left panel of Fig. 4 depicts the standardized linear rank statistics T_µ against tumor size (solid line). The 1st, 10th, 90th, and 99th quantile points of tumor size along with the estimated cutpoint μ̂ of 1.8 are indicated on the axis of tumor size with symbols of black square, black circle, black triangle, and black diamond, respectively. In the top right panel, the CIFs of two groups, tumor size greater than 1.8 (dashed line) and less than or equal to 1.8 (solid line), are displayed. In the same way, plots in the second and third rows of Fig. 4 are respectively drawn with the subgroups of 594 patients having tumor size greater than 1.8 and of 334 patients having tumor size greater than 3. In summary, we could classify the patients with tumor size less than or equal to 1.8 into a low-risk group, while the patients included in a high-risk group (i.e., having tumor size more than 1.8) could be further classified into three subgroups, lowest, moderate, and highest high-risk groups, depending on tumor sizes at surgery of 1.8 to 3, 3 to 6.5, and over 6.5, respectively.