Estimand-based Inference in Presence of Long-Term Survivors
Yi-Cheng Tai Weijing Wang Martin T. Wells

Equivariant Poincaré duality for cyclic groups of prime order and the Nielsen realisation problem
    Kaif Hilmankaif@mpim-bonn.mpg.de Dominik Kirsteinkirstein@mpim-bonn.mpg.de Christian Kremerkremer@mpim-bonn.mpg.de
    September 9, 2024
======================================================================================================================




 



1 Department of Statistics and Data Science, Cornell University, Ithaca, NY, USA

2 Institute of Statistics, National Yang Ming Chiao Tung University, Taiwan, ROC





§ ABSTRACT 
In this article, we develop nonparametric inference methods for comparing survival data across two samples, which are beneficial for clinical trials of novel cancer therapies where long-term survival is a critical outcome. These therapies, including immunotherapies or other advanced treatments, aim to establish durable effects. They often exhibit distinct survival patterns such as crossing or delayed separation and potentially leveling-off at the tails of survival curves, clearly violating the proportional hazards assumption and rendering the hazard ratio inappropriate for measuring treatment effects. The proposed methodology utilizes the mixture cure framework to separately analyze the cure rates of long-term survivors and the survival functions of susceptible individuals. We evaluate a nonparametric estimator for the susceptible survival function in the one-sample setting. Under sufficient follow-up, it is expressed as a location-scale-shift variant of the Kaplan-Meier (KM) estimator. It retains several desirable features of the KM estimator, including inverse-probability-censoring weighting, product-limit estimation, self-consistency, and nonparametric efficiency. In scenarios of insufficient follow-up, it can easily be adapted by incorporating a suitable cure rate estimator.
In the two-sample setting, besides using the difference in cure rates to measure the long-term effect, we propose a graphical estimand to compare the relative treatment effects on susceptible subgroups. This process, inspired by Kendall's tau, compares the order of survival times among susceptible individuals. The proposed methods' large-sample properties are derived for further inference, and the finite-sample properties are examined through extensive simulation studies. The proposed methodology is applied to analyze the digitized data from the CheckMate 067 immunotherapy clinical trial.



Keywords: Cure mixture model, Estimand, Immunotherapy, Insufficient follow-up, Kaplan-Meier estimator, Kendall's tau, Nonproportional hazards, Self-consistency, Sufficient follow-up, Two-sample comparison

Equivariant Poincaré duality for cyclic groups of prime order and the Nielsen realisation problem
    Kaif Hilmankaif@mpim-bonn.mpg.de Dominik Kirsteinkirstein@mpim-bonn.mpg.de Christian Kremerkremer@mpim-bonn.mpg.de
    September 9, 2024
======================================================================================================================





§ INTRODUCTION


In clinical trials for novel cancer therapies, these treatments often exhibit distinct effect profiles compared to conventional treatments such as chemotherapy, likely due to their different pharmacological mechanisms <cit.>. For instance, immunotherapy stimulates the immune system to attack cancer cells, while targeted therapies disrupt cancer cell growth by focusing on specific genetic alterations or proteins, thus minimizing harm to healthy cells. Chemotherapy, often used as the control group in these trials, primarily aims to directly destroy rapidly dividing cancer cells. These varying mechanisms of action lead to different temporal effects on patient outcomes, which are evident in the survival curves. The crossing or delayed separation and subsequent leveling-off in the tails of survival curves, as shown in many reports, clearly violate the proportional hazards assumption.
In such scenarios, the log-rank test may lose power, and the hazard ratio becomes an inappropriate measure for evaluating treatment effects <cit.>.
While novel therapies aim to improve long-term survival, the presence of both long-term survivors and patients who do not respond favorably to treatments introduces complexity into clinical trial analyses <cit.>. The heterogeneity in patient responses often makes interpreting trial results challenging <cit.>.
In accordance with the ICH guidelines <cit.>, formulating suitable estimands for clinical trials involving long-term survivors is crucial. This is essential for the accurate analysis and interpretation of trial results, particularly when traditional statistical methods are inadequate.
 
Empirical evidence indicates that the long-term benefits of promising novel therapies might be limited to a subset of patients <cit.>. This observation aligns with the mixture formulation approach commonly used in analyzing survival data with a cure fraction <cit.>.  
It assumes that there are two groups, one with the potential to be “cured,”  and the other “uncured” or susceptible subjects. In this context, the term 'cure' does not necessarily imply the complete eradication of the disease but may instead represent a long-term survival probability that indicates no further risk of the event of interest. 
Therefore, we adopt this framework to develop suitable estimands for separately measuring the treatment effects in the two sub-populations.

Nonparametric inference of the mixture cure model requires that the follow-up period is sufficiently long to ensure the potential observation of maximum lifetimes in the susceptible sub-population, necessary for the purpose of identifiability <cit.>. In the medical field, landmark trials that aim to establish the long-term efficacy of treatments often include results from extended follow-up periods.
An example is the CheckMate 067 trial, which enrolled 1,296 patients across 137 sites globally and extended its follow-up to at least 6.5 years, providing crucial insights into the long-term outcomes and safety of the treatments <cit.>. On the other hand, statistical methods for scenarios of insufficient follow-up have also been developed. Escobar-Bach and colleagues have utilized extrapolation techniques from extreme value theory to reduce the bias due to missing tail information <cit.>. Furthermore, the paper by Maller et al. reviews recent progress in cure mixture models <cit.>.

In the presence of covariates, various approaches to analyzing the cure mixture model have been explored. Statistical Methods in Medical Research featured a special issue on this topic, edited by Balakrishnan in 2017 <cit.>.
In a two-sample setting, we assess treatment effects by examining differences across key measures. The difference in cure fractions highlights the impact on long-term survivors, while comparisons of survival times among susceptible sub-groups allow us to evaluate the treatment's efficacy for those who do not achieve cure status.
For the latter comparison, the Cox-TEL method, which utilizes a Taylor expansion technique to bridge Cox proportional hazards (PH) and PH cure models for data with long-term survival, adjusts the hazard ratio to more accurately quantify the treatment effect among the two susceptible groups <cit.>. However, the Cox-TEL approach relies on the proportional hazards assumption for the two susceptible groups, an assumption that may be applicable only in certain limited clinical trial scenarios.




In Section <ref>, we first review the cure mixture framework and related existing results. We then focus on a location-scale-shift variant of the Kaplan-Meier (KM) estimator to estimate the susceptible survival function. With sufficient follow-up, we demonstrate that this estimator possesses several useful representations that parallel those of the KM estimator. These include the product-limit form, the Inverse Probability of Censoring Weighting (IPCW) expression, and a self-consistency equation. We also establish the theoretical properties of this latency survival estimator, including weak convergence and nonparametric efficiency. If follow-up is not sufficient, we demonstrate that the susceptible survival function can be estimated by incorporating the cure rate estimator proposed by Escobar-Bach and Van Keilegom <cit.>. Their method involves adding a compensating component to the tail of the KM estimator.

Section 3 examines the cure mixture formulation within a two-sample framework. Beyond assessing long-term survivor outcomes by comparing differences in cure rates between treatment groups, we introduce a graphical estimand that captures the temporal impact of treatments on susceptible groups. The proposed method is a modification of the approach originally developed by Tai et al. <cit.>. This plot offers a clear interpretation of the relative treatment effect on the susceptible groups over time, which is preferred to a single number summary such as an average hazard rate or restricted mean survival rate <cit.>. Our approach does not make any additional assumptions about the relationship between the groups under comparison, enhancing its versatility compared to the Cox-TEL approach.
An extensive simulation study examining the finite-sample properties of the proposed methodology is presented in Section 4, and an application of this methodology to the CheckMate 067 clinical trial is detailed in Section 5. The technical details, proofs of theorems, and additional numerical results are provided in the Supplementary Material.





§ ESTIMATION OF SUSCEPTIBLE SURVIVAL FUNCTION

In the one-sample setting, let T be the failure time with the survival function S(t) = (T > t). Define ξ as the  indicator of susceptibility such that T < ∞ if ξ = 1 and T = ∞ if ξ = 0. Under the mixture formulation, S(t) can be written as 

    S(t) = S_a(t) (1-η) + η,

where  η = Pr(ξ = 0) represents the cure fraction, and S_a(t) = Pr(T > t | ξ = 1) denotes the survival function for susceptible individuals, also known as the latency survival function. Under right censoring, observations that are temporarily censored may become mixed with cured ones (long-term survivors), which can impact the identifiability of η and S_a(t). Let C be the censoring time with the survival function G(t) = (C > t). 
Assume that C is proper, with G(∞) = 0, and that T and C are independent and do not experience simultaneous jumps. Observed variables include X = min(T, C) and δ = I(T ≤ C), where  I(A) is the indicator function that equals 1 if the event A is true and 0 otherwise.  Denote ζ_C and ζ as the right end points of the supports of C and T|ξ = 1, respectively. The condition of sufficient follow-up is related to the requirement that ζ≤ζ_C, which means that the duration of follow-up is long enough to observe the largest event time  among susceptible individuals <cit.>.



 §.§ Nonparametric Analysis of the Cure Mixture Framework: A Review

Denote (T_i,C_i,ξ_i)  (i = 1,…,n) as identically and independently  distributed replications of (T,C,ξ).
Observed variables include (X_i, δ_i), where X_i = T_i ∧ C_i and δ_i = I(T_i ≤ C_i) (i = 1, …, n). Let 0 < t_(1) < … < t_(K) be the  distinct ordered failure times, t_(0) = 0 and K be the number of distinct failure points.  
The KM estimator of S(t) can be written as the product-limit form:

    Ŝ(t) = ∏_k: t_(k)≤ t(1 - d_(k)/y_(k)),

where d_(k) = ∑_i=1^n I(X_i = t_(k), δ_i = 1) and y_(k) = ∑_i=1^n I(X_i ≥ t_(k)) (k = 1, …, K). 
The KM estimator possesses several desirable properties. 
It has been demonstrated that F̂(t) = 1- Ŝ(t) can also be expressed as an IPCW estimator <cit.> such that 

    F̂(t) = ∑_i=1^n I(X_i ≤ t, δ_i = 1)/nĜ(X_i) = 1/n∑_k: t_(k)≤ td̃_(k),

where  d̃_(k) = d_(k)/Ĝ(t_(k)) and Ĝ(t) is the KM estimator of G(t) given by

    Ĝ(t) = ∏_k:u ≤ t(1 - ∑_i=1^n I(X_i = u, δ_i = 0)/∑_i=1^n I(X_i ≥ u)).
 
The mass assigned to t_(k) is given by

    Ŝ(Δ t_(k)) ≡Ŝ(t_(k-1)) - Ŝ(t_(k)) = d_(k)/nĜ(t_(k)) =  1/nd̃_(k),

where Ŝ(t_(0)) = Ŝ(0) =1. 
When the largest observation is a censored observation, Ŝ(t_(K)) remains greater than 0, indicating that the KM curve reaches a plateau or levels off.  

Under sufficient follow-up with ζ≤ζ_C,  the cure fraction η can be estimated by

    η̂ = Ŝ(t_(K))= 1-  ∑_k=1^K d̃_(k)/n.
  
Maller and Zhou established the properties of η̂, such as consistency and asymptotic normality, given specific regularity conditions <cit.>. 
Under insufficient follow-up, Escobar-Bach and Van Keilegom <cit.> propose a new estimator under the assumption that T|ξ = 1 belongs to the maximum domain of attraction of an extreme value
 distribution. To estimate cure rates, they suggest utilizing the tail of the KM estimator, enhanced with a compensating term derived from extrapolation techniques in extreme value theory. Their formula is given by:

    η̌_b = η̂ - Ŝ(b t_(K)) - Ŝ(t_(K))/b̌_γ - 1,

where b ∈ (0,1) is a scaling factor,  γ is an extreme value index that controls the tail behavior of S_a(t) and

    b̌_γ = Ŝ(b t_(K)) - Ŝ(b^2 t_(K))/Ŝ(t_(K)) - Ŝ(b t_(K)).

Consequently, the survival function can be estimated by directly substituting an estimator of η, either η̂ or η̌_b, into equation (<ref>). The referenced paper suggests using bootstrap resampling to select the value of b <cit.>. 








 §.§ The Susceptible Survival Estimator under Sufficient Follow-Up

Under sufficient follow-up, substituting η̂ for η in equation (<ref>) and applying it to equation (<ref>) yields the estimator:

    Ŝ^LS_a(t) = Ŝ(t) - η̂/1- η̂,

which can be interpreted as a location-scale-shift variant of Ŝ(t).
We now demonstrate that Ŝ^LS_a(t) in (<ref>) is equivalent to the alternative IPCW and product-limit estimators:

    Ŝ^W_a(t)= ∑_i=1^n I(X_i > t, δ_i = 1)/n(1-η̂)Ĝ(X_i) = 
      1/n̂_a∑_k: t_(k) >td̃_(k)
      = ỹ_(k)/n̂_a

and

    Ŝ_a^PL (t) = ∏_k: t_(k)≤ t(1 - d̃_(k)/ỹ_(k)),

where n̂_a = n(1-η̂), and ỹ_(k) is the adjusted number at risk at t_(k) given by

    ỹ_(k)=∑_i=1^n I(X_i ≥ t_(k),δ_i=1)/Ĝ(X_i-) = ∑_j=k^Kd̃_(k).

Note that ỹ_(1) = ∑_k=1^Kd̃_(k)= n̂_a.  The expression Ŝ_a^W(t) in (<ref>) is an  IPCW  estimator of S_a(t), with the mass assigned to each observed failure point adjusted by the estimated sample size n̂_a for the susceptible group. In addition, Ŝ_a^PL (t) in (<ref>) is an adjusted variant of the product-limit estimator in (<ref>).

Based on mass calculations, it can be shown that the three expressions in (<ref>)-(<ref>) are identical. Specifically, from the product-limit expression in (<ref>), the mass at t_(k+1) is given by:

    Ŝ_a^PL(Δ t_(k+1)) = Ŝ_a^PL (t_(k))  d̃_(k)/ỹ_(k).

Since Ŝ^PL_a(t_(0)) =Ŝ^PL_a(0) =1, Ŝ^PL_a(Δ t_(1)) = d̃_(1)/ỹ_(1). It is straightforward to show that Ŝ^PL_a(t), Ŝ^W_a(t) and Ŝ^LS_a(t) are equal and can be written as

    Ŝ_a(t) = ∑_i=1^n I(X_i > t, δ_i = 1)/n(1-η̂)Ĝ(X_i) 
    = 1/n̂_a∑_k: t_(k)>td̃_(k) = ỹ_(k)/n̂_a.

We observe that the effect of long-term survivors on the estimation of S_a(t) simply involves reducing the sample size to the susceptible subgroup. Consider the following three risk sets:

    R(t)    = {i: X_i ≥ t | i=1,…,n}, 
    
    R_a(t)    = {i: X_i ≥ t, ξ_i=1 | i=1,…,n}, 
    R(t)    = {i: X_i ≥ t, δ_i=1 | i=1,…,n}.

When there are long-term survivors with ξ_i=0, not all members in R(t) are susceptible to the event of interest. However, in the presence of censoring, the susceptible risk set R_a(t) is not fully recoverable from the available data. Given that R(t) ⊂ R_a(t) ⊂ R(t), R(t) can be used as a proxy for R_a(t). To account for selection bias among observations in R(t), the inverse-probability-of-censoring weighting technique is employed to estimate the hazard rate of T|ξ = 1 at time t_(k) as d_(k)/y_(k), which explains (<ref>). 


The theoretical properties of Ŝ_a(t)  outlined in this subsection are instrumental for extending the analysis to different censoring scenarios, such as interval censoring, and offer a rationale for employing bootstrap re-sampling techniques for additional inferential objectives.
In a remarkable article, Gill <cit.> showed that the nonparametric maximum likelihood estimator  is often determined as the solution of the likelihood equations for a collection of smooth parametric submodels. These equations are, in fact, precisely the “self-consistency” equations introduced by Efron <cit.>. Consequently, in many settings, a solution of the self-consistency equation is equivalent to the nonparametric maximum likelihood estimator, and in general settings, the nonparametric maximum likelihood estimator will be asymptotically efficient. 
Following the approach of Strawderman and Baer <cit.>, we derive the self-consistency property for 
Ŝ^PL_a(t) in the cure model setting. The proposed self-consistency equation is a data-dependent function that represents the expected number of observations which are both susceptible and surviving at each time point in the sample. This function is expressed in terms of the target function S_a(t)  and specific nonparametric estimates. The representation of the self-consistency function, along with its proof, is provided in Section 1 of the Supplementary Materials. The adaptability of the self-consistency property to interval-censored data and more complex censoring schemes presents a potential avenue for future extensions within the cure mixture framework
<cit.>.
















































We further study the weak convergence and nonparametric efficiency of the susceptible survival function. The weak convergence property of Ŝ_a(t) is given in  Theorem <ref> and the property of nonparametrically efficiency is stated in  Theorem <ref>. 

    Assume that 0 ≤η < 1,
    
    sup_0≤ t < ∞|Y(t)/n - S(t)G(t)|=o_P(1)
 and
    
    ∫_0^ζdF(u)/G(u-)<∞.

    Under sufficient follow-up with ζ≤ζ_C, √(n)(Ŝ_a(t) - S_a(t)) converges weakly to a mean-zero Gaussian process with t∈[0,ζ].
    

Note that when ζ < ζ_C, the condition in equation (<ref>) is satisfied. Therefore, it is necessary to verify condition (<ref>) only when ζ=ζ_C.  The sufficiency of the condition in (<ref>) is discussed in detail by Gill <cit.>.



The theory of regular and asymptotically linear (RAL) estimators can be applied to analyze Ŝ_a(t). This estimator is represented in terms of the efficient influence function, a concept introduced in the referenced book <cit.>. 
The influence function of the RAL estimator, which possesses the lowest asymptotic variance (van der Vaart, 1998, Theorem 25.20), corresponds to that of the asymptotically efficient estimator <cit.>.



Recall that the estimators Ŝ^PL_a(t), Ŝ^W_a(t), and Ŝ^LS_a(t) all equal Ŝ_a(t) in (<ref>).
The efficient influence function form of a survival function estimate for general right-censored data is detailed in Section 4 of the book by Robins and Rotnitzky <cit.> and on page 133 of the book by Van der Laan and Robins <cit.>.
Specializing those results for the model in (<ref>) and the estimation of η gives the efficient influence function of Ŝ_a(t). 











    ψ_t(X,δ)=I(X>t,δ=1)/(1-η)G(Y-)+∫_0^ζq_t(u)/(1-η)π(u)dM_C,1(u)-S_a(t),

where π(u)=lim_n→∞ Y(u)/n, 

    q_t(u) = lim_n →∞∑_i=1^n I(X_i>t,δ_i=1,X_i≥ u)/{nG(X_i-) },

 and 

    M_C,i(t)=I(X_i≤ t,δ_i=0)-∫_0^tI(Y(u)>0)-dG(u)/G(u-).










 

    Under sufficient follow-up with ζ≤ζ_C, Ŝ_a(t) is regular and asymptotically linear (RAL) with the influence function

    ψ_t(X, δ) = δ{I(X>t) - S_a(t)}/(1-η)G(X)  + ∫_0^ζq_t(u)/π(u) dM_C,1(u),

where π(u) and q_t(u) are the probability limits of Y(u)/n and 

    ∑_i=1^nδ_iI(X_i≥ u) [I(X_i>t) - S_a(t)]/n(1-η)G(X_i),

respectively, λ_c(t) is the hazard function of C, and the censoring martingale is

    M_C,i(t) =  I(X_i ≤ t, δ_i = 0) - ∫_0^t I(X_i ≥ u) λ_c(u)du.

It follows that Ŝ_a(t) is nonparametrically efficient and √(n)(Ŝ_a(t) - S_a(t)) converges in distribution pointwise to a mean-zero normal random variable for t∈[0,ζ].



The representation in Theorem <ref> is reminiscent of the representation of the KM estimator process as an identically and independently distributed process, as given by Lo and Singh <cit.>. This representation justifies the bootstrap method for estimating the standard error of functionals of the  susceptible survival function estimate and its quantiles. Based on the bootstrap resampling approach, it provides a way of constructing confidence intervals (bands) for the unknown parameters (functionals of the distribution or quantile function). The proofs of the two theorems above are provided in Subsections 2.2 and 2.3 of the Supplementary Material. 



 §.§ The Susceptible Survival Estimator under insufficient Follow-Up

Under conditions of insufficient follow-up, by substituting η̌_b from equation (<ref>) into equation (<ref>), we obtain:

    Š_a^LS(t;b) = Ŝ(t) - η̌_b/1 - η̌_b.

The properties of Š_a^LS(t;b) are determined by those of Ŝ(t) and η̌_b. For detailed properties of η̌_b, refer to Escobar-Bach and Van Keilegom <cit.>. Although η̌_b does not ensure convergence to the true cure rate η, it helps mitigate the underestimation typically associated with η̂ under insufficient follow-up.
As the duration of follow-up increases and approaches ζ, the estimator η̌_b becomes more accurate in approximating the true cure rate η. This improvement in the accuracy of η̌_b directly enhances the performance of Š_a^LS(t;b), making it more reliable in estimating the survival function S_a(t). In the supplementary materials, we offer heuristic discussions on the properties of Š_a^LS(t;b), with a particular focus on analyzing its bias relative to S_a(t).




§ TWO-SAMPLE APPLICATION
 


We propose a flexible alternative by adapting the methodology of Tai et al. <cit.>, which does not rely on specific model assumptions such as proportional hazards. Specifically, the tau process, which measures the relative performance of two groups under comparison, is defined as follows:

    τ(t) = ∫_0^t S_1(u) dF_0(u) - ∫_0^t S_0(u) dF_1(u).

Each component of the integrand of τ(t) represents the process where, at each failure time in one group, the survival probability of the other group is evaluated at that specific time point. The function τ(t) sums these differences in survival probabilities up to time t, providing a cumulative measure of the disparity between the two groups over time.  A positive value of τ(t) indicates that the treatment group (Group 1) exhibits a better effect up to time t. Note that τ(t) is a unitless and model-free measure, which makes it a robust and interpretable treatment-effect estimand, especially in the presence of non-proportional hazards <cit.>. Let τ̂(t) denote the estimator of τ(t). 

In the absence of a cure, such that (T < ∞) = 1, τ(∞) represents Kendall's tau correlation between the group indicator and the failure time. When cure is a possibility, the cross-comparison via the difference in the integrands of τ(t) is blurred by the presence of long-term survivors, as S_1(u) and S_0(u) may plateau at different levels. Specifically, when Group 1 exhibits a much higher cure rate than Group 0, leading to S_1(u) ≫ S_0(u) at large values of u, the differences in early stages may become obscured. In the next subsection, we define another tau process to compare event times between the two susceptible subgroups. We then explore its relationship with τ(t), as well as the cure rates η_0 and η_1.





 §.§ Tau Process for Susceptible Subgroups



Under the two-sample setting, let (T_ℓ,C_ℓ,ξ_ℓ) be the failure time, censoring time and the indicator of susceptibility for Group ℓ with ℓ=0,1. The mixture model described in (<ref>) is adapted for each subgroup with ξ_ℓ = 1, such that

    S_ℓ(t)  = S_a,ℓ(t) (1-η_ℓ) + η_ℓ,

where  S_ℓ(t) = (T_ℓ > t),  S_a,ℓ(t) = (T_ℓ > t | ξ_ℓ = 1)  and η_ℓ = (ξ_ℓ = 0) for ℓ = 0,1. Define F_a,ℓ(t) = 1- S_a,ℓ(t) as the distribution function of T_j|ξ_j = 1 (j = 0,1). The long-term treatment effect can be described by the difference of η_1 and η_0. The treatment effect on the susceptible groups is evaluated by quantifying the difference between S_a,1(t) and S_a,0(t). 

To characterize the treatment effect on the two susceptible groups, we introduce the susceptible tau process as follows:

    τ_a(t)=∫_0^tS_a,1(u)dF_a,0(u)-∫_0^tS_a,0(u)dF_a,1(u).

For susceptible patients who do not ultimately achieve a cure, a positive value of τ_a(t) suggests that the treatment may still prolong the time until the occurrence of the unfavorable event of interest. Given that τ(t) = (T_0 < T_1 ∧ t) - (T_1 < T_0 ∧ t), 
we may think of τ_a(t) as follows:

    τ_a(t) = 1/(1-η_0)(1-η_1)[ ((T_0 < T_1 ∧ t) - η_1 F_0(t)) - ((T_1 < T_0 ∧ t) - η_0 F_1(t)) ],

where -η_1 F_0(t) reflects the adjusted probability, excluding comparisons between susceptible subjects who die before t in group 0 and long-term survivors in group 1;  -η_0 F_1(t) reflects the adjusted probability, excluding comparisons between susceptible subjects who die before t in group 1 and long-term survivors in group 0; and {(1-η_0)(1-η_1)}^-1 serves as a normalizing constant.


Note that the signs of τ(t) and τ_a(t) may differ, indicating a potential reversal in effects between susceptible individuals and those considered long-term survivors. Specifically, we can write 

    τ(t)=(1-η_0)(1-η_1)τ_a(t)+(1-η_0)η_1F_a,0(t)-(1-η_1)η_0F_a,1(t).

Consider a situation that η_1 - η_0 is far greater than zero but S_a,1(t) < S_a,0(t) (equivalently, F_a,1(t) > F_a,0(t)). Equation (<ref>) indicates that it is still possible for τ(t) > 0 while τ_a(t) < 0. In certain oncology clinical trials, this scenario may correspond to a situation where a subset of patients in Group 1 does not respond favorably to immunotherapy, resulting in their classification as short-term survivors, with even shorter survival. 

It has been argued that high-risk individuals are more likely to be depleted over time, which can make the hazard ratio susceptible to selection bias <cit.>. Additionally, differences in cure rates may not fully capture treatment effects in the early stages of a trial, which are crucial for assessing overall efficacy. Therefore, to provide a comprehensive assessment of treatment effects, we recommend using multiple measures, including τ(t), τ_a(t), and η_1 - η_0.




 §.§ Estimation of τ_a(t)

Observed data in the two-sample setting can be denoted as (X_ℓ,i, δ_ℓ,i) for  i = 1, …, n_ℓ and ℓ=0,1.
Let X̃_ij = X_0,i∧  X_1,j and O_ij=I(X_0,i<X_1,j,δ_0,i=1)+I(X_0,i>X_1,j,δ_1,j=1). Under sufficient follow-up, η̂_ℓ is a legitimate estimator of η_ℓ for ℓ = 0,1.  The resulting estimator of τ_a(t) is given by 

    τ̂_a(t) = ∑_i,jψ̂_ij(t)ŵ_ij/n_0 n_1(1-η̂_0)(1-η̂_1),

where

    ψ̂_ij(t) = Õ_ijsign(X_1,j-X_0,i)I(X_ij≤ t)/Ĝ_0(X_ij)Ĝ_1(X_ij),

and

    ŵ_ij = [(1 - η̂_0)Ŝ_a,0(X_0,i)/(1 - η̂_0)Ŝ_a,0(X_0,i) + η̂_0]^(1 - δ_0,i)×[(1 - η̂_1)Ŝ_a,1(X_1,j)/(1 - η̂_1)Ŝ_a,1(X_1,j) + η̂_1]^(1 - δ_1,j).

The asymptotic property of τ̂_a(t) is stated in Theorem <ref>, and its asymptotic variance can be obtained using the bootstrap approach, with the proof provided in Subsection 2.3 of the Supplementary Material.
 


    Suppose that the conditions of Theorem <ref> hold for both groups, and let n_1/n converge to p_1 as n_0 and n_1 tend to infinity, where 0 < p_1 < 1. If E[ψ_ij(t)w_ij]^2 is finite for t up to min(ζ_0, ζ_1), then √(n)(τ̂_a(t) - τ_a(t)) converges pointwise to a mean-zero normal random variable.

One can apply a bootstrap approach for right censored data to develop inferential procedures.   The validity of such bootstrap confidence interval and tests follow using classical arguments for right censored data <cit.>. The processes τ̂(t) and τ̂_a(t) can be implemented using the R package 'tauProcess' <cit.>.

Under insufficient follow-up, the estimator of η_ℓ can be obtained by modifying Equation (<ref>) for the two-sample setting, and is denoted as η̌_ℓ, b_ℓ for ℓ = 0, 1.
The corresponding estimator of τ_a(t) can be modified by replacing η̂_ℓ with η̌_ℓ, b_ℓ, and is denoted as τ̌_a(t; b_0, b_1).
Note that the tail behaviors of the two groups may differ, necessitating the separate estimation of b_0 and b_1 in practical applications. The implementation procedure suggested by Escobar-Bach and Van Keilegom will be summarized in the data analysis section <cit.>.








§ SIMULATION STUDY
 
In the first design, we assess the performance of the estimators for S_a(t) and η under both sufficient and insufficient follow-up conditions. In the second simulation design, we evaluate the estimators of τ_a(t) for two scenarios—crossing and non-crossing survival functions S_a,0(t) and S_a,1(t), with η_ℓ values of 0.2 and 0.4 for ℓ = 0, 1. For both simulation settings, we analyze the average bias of each estimator, along with the bootstrap standard deviation estimates and the empirical coverage probabilities of the estimated confidence intervals. These results are derived using 2000 bootstrap resamples across 500 simulation runs. Additional simulation results are detailed in the Supplementary Material.





 §.§ Finite-Sample Performance of Estimators for S_a(t) and η


In the initial design, we assess the finite-sample performance of estimators for S_a(t) and η, setting η = 0.2 and 0.4. To model T|ξ = 1 with bounded support, we employ a Beta distribution with parameters α_1 = 1 and α_2 = 3.
We first discuss the results under conditions where ζ = ζ_C and (<ref>) is met, with the censoring variable C following a Uniform[0, 4] distribution. As shown in Table <ref>, the average biases of Ŝ_a(t) and η̂ are close to zero, and the bootstrap estimates of the standard deviation for Ŝ_a(t) closely match the empirical estimates. 
The empirical coverage probabilities are also around the 95% nominal level. Note that confidence intervals are wider under η = 0.4 compared to η = 0.2. Additional analyses based on other scenarios are presented in the Supplementary materials.



[Insert Table <ref> & Table <ref>]




Under insufficient follow-up,  the censoring variable C following a Uniform[0, 0.8] distribution.
We evaluate the performance of Š_a^LS(t;b_*) using η̌_b_* to estimate η, where b_* is selected such that the corresponding estimator η̌_b closely matches the average outcome of a bootstrap experiment <cit.>. This choice minimizes the deviation between η̌_b and the average bootstrap estimator across multiple resampled datasets.
Compared with the performance of η̂ shown in Table <ref>, Š_a^LS(t;b_*) exhibits slight bias, as indicated in row (a) of the last column in Table <ref>.
Nevertheless, the performance of Š_a^LS(t;b_*) remains satisfactory as an estimator of S_a(t), although the coverage probability tends to deviate more from 95% as t increases.  





 §.§ Finite-Sample Performance of τ̂_a(t)

In this subsection, we present results for τ̂_a(t) under conditions of sufficient follow-up. Note that under insufficient follow-up, the performance of τ̌_a(t; b_0, b_1) is influenced by the shapes of the unobserved tail distributions. Therefore, we apply this modified estimator in our data analysis but exclude it from the simulations.
Figures S4 and S5 in the Supplementary material depict  two non-crossing susceptible survival functions with the same cure rates and the corresponding τ_a(t). We present the results for the case with apparent disparity, and the other case is provided in the Supplementary material.
From Table <ref>,  
we observe that the average bias of τ̂_a(t) is almost zero. The bootstrap estimates of the standard deviation of τ̂_a(t) closely align with the empirical estimates, and the empirical coverage probabilities are close to the 95% nominal level.  
Notice that the lengths of the confidence intervals are associated with the values of η_0 and η_1. Higher cure rates lead to wider confidence intervals for τ_a(t). According to Table <ref>, τ̂_a(1) yields a significant result for testing H_0: τ_a(1) = 0.



[Insert Table <ref>]



In the Supplementary Material, Figures S6 and S7 depict scenarios where the susceptible survival functions intersect at a susceptible survival probability of approximately 0.5. The former corresponds to η_0=η_1 = 0.2, and the latter corresponds to η_0=0.2, η_1 = 0.4. The corresponding simulation results are summarized in Table <ref>. The performance of τ̂_a(t) is similar to that of the previous settings. The result from the last column indicates that there is no significant result for testing H_0: τ_a(1) = 0 based on τ̂_a(1).



[Insert Table <ref>]



In Figure S8 of the Supplementary Material, we present simulation results that examine the situation where η_1 > 0 and η_0 = 0. These results confirm the validity of the inference procedure based on τ̂_a(t). 





§ DATA ANALYSIS
 
The CheckMate 067 trial, a randomized, multicenter, phase 3 study conducted from July 2013 to March 2014, evaluated three different immunotherapy strategies for advanced melanoma in a cohort of 945 patients. Participants were assigned to receive either a combination of nivolumab and ipilimumab, nivolumab alone, or ipilimumab alone, with each regimen designed to boost the immune system's ability to combat the disease.

First, we apply our proposed methodology by digitizing the overall survival KM curves for the three treatment groups based on the 4-year report <cit.>.
The estimated cure rates based on the tail values of the KM curves are 0.516 for the combination treatment, 0.452 for nivolumab alone, and 0.272 for ipilimumab alone. Given that the 4-year data may not represent sufficient follow-up, we also applied the extrapolation method proposed by Escobar-Bach and Van Keilegom <cit.>. This approach provided cure rate estimates based on η̌_b_*, yielding values of 0.481 for combination therapy, 0.421 for nivolumab alone, and 0.199 for ipilimumab alone. 
The difference in cure rates between the combined group and 'ipilimumab alone' is significant, calculated as 0.516 - 0.272 = 0.244 (0.164, 0.323) with p-value 1.93× 10^-9 using the KM tail estimates, and 0.481 - 0.199 = 0.282 (0.128, 0.445) with p-value 3.97× 10^-4 using the modified estimator. Figure <ref> offers a graphical assessment, showcasing the KM curves, τ̂(t), as well as both versions of the estimated susceptible survival functions and the susceptible tau process. In plot (A), the combined group shows a higher survival curve and exhibits a higher cure rate. Plot (B) shows that τ̂(t) becomes positive starting from the third month onward.
Using η̂_ℓ as the estimate of η_ℓ, plots (C) and (D) display Ŝ_a,ℓ(t) for ℓ = 0, 1 and τ̂_a(t), respectively.
Using η̌_ℓ,b_*,ℓ as the estimate of η_ℓ, plots (E) and (F) display Š_a,ℓ^LS(t;b_*,ℓ) for ℓ = 0, 1 and τ̌_a(t; b_*,0, b_*,1), respectively. By comparing the last two rows of Figure <ref>, we observe that the modified cure rate estimators do not significantly impact the estimation of the susceptible survival functions and the susceptible tau process.
Although the combined therapy clearly outperforms 'ipilimumab alone' based on the KM estimators and τ̂(t), the analysis of susceptible groups indicates  reversal relationships. This suggests that while a larger proportion of patients in the combined therapy group are long-term survivors, those who do not achieve a durable effect have similar or slightly shorter survival times compared to their counterparts in Group 0.
The difference in cure rates between 'nivolumab alone' and 'ipilimumab alone' is significant, calculated as 0.452 - 0.272 = 0.18 (0.101, 0.260) with p-value 9.08 × 10^-6 using the KM tail estimates, and 0.421 - 0.199 = 0.222 (0.067, 0.384) with p-value 0.005 using the modified estimator. According to Figure <ref>, although the 'nivolumab alone' group exhibits a higher cure rate, the estimated susceptible survival functions are roughly similar. The estimated curves of τ_a(t) demonstrate their proximity to zero and even suggest a slight reversal, as observed in the third row.





[Insert Figure <ref> & Figure <ref>]
















§ CONCLUDING REMARKS
 
Long-term survivorship, often referred to as cure, is a common outcome in various fields. In the context of cancer treatment, the sustained effectiveness of immunotherapy and other advanced options has generated optimism among both patients and healthcare professionals. 
However, understanding the potential heterogeneity in patient responses to these treatments remains a crucial aspect of devising appropriate treatment strategies for the right individuals. Cure mixture models provide a useful framework that allows for separate evaluation of long-term survivors and susceptible individuals who do not achieve the desired long-term status.
While the estimation of the cure rate assumes sufficient follow-up, which may not always be attainable, extreme value theory provides methods to develop new estimators <cit.>. These estimators can reduce bias in the tail estimates of the KM curves under conditions of insufficient follow-up.
Under sufficient follow-up, we demonstrate that the estimator of the susceptible survival function, Ŝ_a(t), inherits many favorable properties of the KM estimator, which paves the way for further extensions across various data structures and settings.
Under insufficient follow-up, the location-scale-shift version adapts the modified cure rate estimator suggested by Escobar-Bach and Van Keilegom <cit.>. To evaluate the effect on long-term survivors, we recommend using the difference in cure rates. Additionally, our proposed graphical estimand, τ_a(t), offers insights into the treatment effects over time for individuals who have not been cured, highlighting the timing and impact of the therapy. The susceptible tau process can be estimated nonparametrically, provided that a suitable cure rate estimator is available.
Building on these methods, employing multiple estimands in clinical trial analysis provides a clear framework that enhances understanding of survival outcomes. This approach not only enriches the analysis but also informs future therapeutic strategies and research directions. 



§ ACKNOWLEDGEMENTS

We express our gratitude to Robert Strawderman and Benjamin Baer for their invaluable discussions on self-consistency and efficiency theory. Additionally, we are thankful to Mikael Escobar-Bach and Ingrid Van Keilegom for kindly providing the code for their cure rate estimator.



    
The authors thank Robert Strawderman and Benjamin Baer for their helpful discussions on self-consistency and efficiency theory. 



The authors declare no potential conflict of interests.




Wang's research received support from the National Science and Technology Council of Taiwan through grants 111-2118-M-A49-008 and 112-2118-M-A49-002. Wells’ research was partially supported by National Institutes of Health awards and R01GM135926 1P01-AI159402.
















    











    































unsrt



