Nonparametric Analysis of Delayed Treatment Effects using Single-Crossing Constraints
aa r X i v : . [ s t a t . M E ] J a n Nonparametric Analysis of Delayed TreatmentEffects using Single-Crossing Constraints
Nicholas C. Henderson ∗ , Kijoeng Nam and Dai Feng Department of Biostatistics, University of Michigan, Ann Arbor, MI, USA BARDS, Merck & Co., Inc., North Wales, PA, USA Data and Statistical Sciences, AbbVie Inc., North Chicago, IL, USA*email: [email protected]
Abstract
Clinical trials involving novel immuno-oncology (IO) therapies frequently exhibit survival pro-files which violate the proportional hazards assumption due to a delay in treatment effect, andin such settings, the survival curves in the two treatment arms may have a crossing before thetwo curves eventually separate. To flexibly model such scenarios, we describe a nonparametricapproach for estimating the treatment arm-specific survival functions which constrains these twosurvival functions to cross at most once without making any additional assumptions about howthe survival curves are related. A main advantage of our approach is that it provides an estimateof a crossing time if such a crossing exists, and moreover, our method generates interpretablemeasures of treatment benefit including crossing-conditional survival probabilities and crossing-conditional estimates of restricted residual mean life. We demonstrate the use and effectivenessof our approach with a large simulation study and an analysis of reconstructed outcomes from arecent combination-therapy trial.
Key words : censored data; clinical trial; constrained estimation; immuno-oncology; non-proportional haz-ards
1. Introduction
Recent advances in immuno-oncology (IO) therapies for the treatment of cancer have led tothe development of a variety of treatments which show great potential for improving long-termpatient outcomes. While very promising, patient response to such immunotherapies is often quitedifferent when compared to more traditional cytotoxic agents such as chemotherapy. Indeed, it iswell-recognized that IO drugs frequently exhibit a clear delay in treatment effect when they arecompared with standard chemotherapies, and moreover, the nature of the delay in this treatmenteffect is often such that the estimated survival curves in the two treatment arms have a crossingat some time point after randomization. Due to this feature of immunotherapies, traditionalcomparisons between IO drugs and chemotherapies can have a number of limitations.For time-to-event endpoints in randomized clinical trials, the log-rank test is the conventionalchoice for testing the superiority or non-inferiority of an active treatment over a control, andwhenever the log-rank test passes a threshold for statistical significance, a hazard ratio is typicallyreported as one of the chief measures of treatment efficacy. However, in settings with survivalcurve crossings where the proportional hazards assumption is plainly violated, the interpretation1
N.C. Henderson & K. Nam & D. Feng of both the log-rank test and an estimated hazard ratio from a Cox proportional hazards modelcould be unclear. A variety of alternatives to the hazard ratio have been suggested in the contextof evaluating the efficacy of immunotherapies or other context where the proportional hazardsassumption regularly fails. These include, for example, difference in restricted mean survival time(RMST) (Zhao and others (2016), Pak and others (2017)), difference in milestone survival (Chen(2015)), average hazard ratios (Schemper and others (2009)), and the proportion of patients thatare functionally “cured” (Chen (2013)). In the context of testing the equality of the arm-specificsurvival curves under possible delayed treatment effects, a number of weighted log-rank tests havebeen suggested including, for example, the piecewise weighted log-rank test (APPLE) proposedin Xu and others (2017) and weighted log-rank tests from the the Fleming-Harrington familywhich can place more weight on later time points (Harrington and Fleming (1982), Rahman andothers (2019)). Combination tests which combine multiple weighted log-rank test statistics havealso been proposed for analyzing IO trials (Lin and others (2020)).While the measures of treatment effect listed above have a clear interpretation in the ab-sence of proportional hazards and reporting such measures of treatment efficacy can certainlybe valuable, imposing estimation constraints on how the survival curves may cross can provideadditional information by which to evaluate a treatment that exhibits delayed benefit. To thisend, we propose nonparametric estimates of the treatment-specific survival functions that allowfor at most one crossing of the two survival functions without requiring that the crossing time bepre-specified. The approach we describe for estimating the survival curves under a single-crossingconstraint can be thought of as a two-stage procedure. In the first stage, one finds conditionalestimates of the survival function by maximizing a nonparametric log-likelihood that is condi-tional on the value of two crossing parameters. Then, in the second stage, one estimates thesecrossing parameters by maximizing a profile log-likelihood function. This produces estimates ofthe two survival curves which satisfy the single-crossing constraint, and it generates an estimateof the crossing time and an estimate of which survival curve is initially dominant. While flexi-ble semiparametric approaches allowing for the crossing of treatment-specific survival functionshave been proposed (see, e.g., Yang and Prentice (2005) and Demarqui and Mayrink (2019)),our approach makes no assumptions about how the treatment-specific survival curves are relatedother than that they cross at most once, and moreover, our approach provides a likelihood-basedestimate of the crossing time if a positive crossing time is determined to be present.One of the advantages of imposing a single-crossing constraint is that situations where thereis a delay in treatment effect and improved long-term survival in the active treatment arm areoften better modeled as having a single crossing after which point the survival curve in the activetreatment arm consistently dominates the survival curve of the control arm. Such patterns ofdelayed treatment effect have been observed in many recent immuno-oncology trials, and hence,constrained modeling in such trials has the potential to improve estimation performance andprovide greater interpretability when comparing the two arm-specific estimates of the survivalfunctions. An important advantage of using single-crossing constraints is that it yields estimatesand uncertainty intervals for the time that the crossing occurs. The estimate of the crossing-timeparameter can be useful in both assessing when the active treatment begins to show superiorityand as an interpretable measure that can be used as a component of additional measures oftreatment efficacy that are designed for scenarios with delayed treatment effect which we outlinein detail in Section 3. These efficacy measures include the proportion of patients surviving up tocrossing, crossing-time conditional restricted residual mean life, crossing-time conditional survivalprobabilities, and pre/post-crossing average hazard ratios.This paper is organized as follows. In Section 2, we describe a two-stage nonparametric esti-mation procedure which assumes that the treatment arm-specific survival curves cross at most elayed Treatment Effects with Single-Crossing Constraints
2. Nonparametric Estimation for Survival Curves with Single Crossings
Survival Curve Profiles and Notation
We assume that n patients have been enrolled in a randomized clinical trial consisting of twotreatment arms. The primary outcome is the time to some event of interest, and for the i th individual in the study, we let T i denote the time-to-failure for this event of interest. Instead ofdirectly observing T i , we observe the follow-up time Y i = min { T i , C i } and the event indicator δ i = I ( T i C i ) where C i denotes the censoring time, and I ( · ) denotes the indicator function.We let A i = 1 denote that patient i was assigned to the active treatment arm, and we let A i = 0 denote that patient i was assigned to the control arm. We also assume that censoringis noninformative in the sense that the terms from the censoring distribution can be factoredout of the likelihood function (Lawless (2011)), and we assume that C i and the treatment armassignment A i are independent.We let S a ( t ) = P ( T i > t | A i = a ) denote the survival function for those patients assigned totreatment group a ∈ { , } . We consider an analysis where these survival functions are allowedto exhibit four distinct “survival profiles” according to whether or not a survival curve crossingoccurs. Specifically, we allow for the possibility that the two survival functions S ( t ) and S ( t )“cross”, but we limit the number of crossings so that they can occur at most once. The crossingrestriction implies that we will have one of four survival profiles, where each profile refers to adistinct crossing pattern of the survival functions. These four possible survival profiles are asfollows: (1) a survival profile where treatment arm a = 1 completely dominates a = 0, namely S ( t ) > S ( t ) for all t >
0; (2) a survival profile where a = 0 dominates a = 1 before somecrossing time θ but a = 1 dominates a = 0 afterwards, i.e., S ( t ) > S ( t ) for 0 t θ and S ( t ) > S ( t ) for t > θ ; (3) a survival profile where a = 1 dominates a = 0 before some crossingtime θ but a = 0 dominates a = 1 afterwards; and (4) a survival profile where a = 0 completelydominates a = 1. Figure 1 illustrates each of these four possible survival profiles.To enforce one of the four possible profiles implied by Figure 1, we introduce the two crossingparameters θ > γ ∈ {− , } . The parameter θ represents the crossing time of the survivalfunctions, and the discrete parameter γ determines which treatment arm has the initially domi-nant survival function. Specifically, when γ = 1 the survival functions are assumed to exhibit thefollowing behavior S ( t ) > S ( t ) for 0 t θ and S ( t ) S ( t ) for t > θ, (2.1)and when γ = −
1, the survival functions are assumed to obey the following inequalities S ( t ) S ( t ) for 0 t θ and S ( t ) > S ( t ) for t > θ . (2.2)If there is a θ > θ is a non-trivial crossing time of the survival functions S and S . If there is no such θ >
0, we
N.C. Henderson & K. Nam & D. Feng S u r v i v a l P r obab ili t y TreatmentControl g = 1 q = 0 Time 0.00.20.40.60.81.0 S u r v i v a l P r obab ili t y TreatmentControl g = 1 q = 2.75 Time0.00.20.40.60.81.0 S u r v i v a l P r obab ili t y TreatmentControl g = −1 q = 0 Time 0.00.20.40.60.81.0 S u r v i v a l P r obab ili t y TreatmentControl g = −1 q = 2.75 Time
Fig. 1. Four possible survival profiles. The crossing parameters θ and γ refer to the crossing-time pa-rameter and initial dominance parameter respectively. We set θ = 0 if no crossing occurs. A value of γ = 1 indicates that the control treatment arm is dominant before the crossing while a value of γ = − must either have S ( t ) > S ( t ) or S ( t ) > S ( t ) for all t > θ to the trivialcrossing time θ = 0.When constructing estimates of the survival functions, we can only enforce the infinite-dimensional constraints implied by (2.1) and (2.2) at a finite number of time points. In ourimplementation, we enforce the survival function constraints (2.1) and (2.2) at each of the ob-served event times. To this end, we let 0 < t < t < ... < t m denote the unique, ordered eventtimes from both treatment arms, and hence, m represents the total number of events occurring ineither of the treatment arms. For γ = 1 and a fixed value of θ , we enforce the following constraints S ( t j ) > S ( t j ) for all j such that t j θS ( t j ) S ( t j ) for all j such that t j > θ . (2.3) elayed Treatment Effects with Single-Crossing Constraints γ = − θ , we enforce the following constraints S ( t j ) S ( t j ) for all j such that t j θS ( t j ) > S ( t j ) for all j such that t j > θ . (2.4)Note that (2.3) and (2.4) together imply that we will have a different set of constraints for eachchoice of ( θ, γ ). In the next subsection, we describe our procedure for estimating the arm-specificsurvival functions when both θ and γ are assumed to be fixed.2.2 Estimating the Survival Functions with Fixed Crossing Parameters
We first describe nonparametric estimation of the survival functions S a ( t ) assuming both θ and γ are known. As detailed in Park and others (2012) in the context of finding the constrainednonparametric maximum likelihood estimate of survival functions under a stochastic orderingconstraint, the constrained maximum likelihood estimates of the survival functions will be discretewith potential jumps only at the event times t , . . . , t m . Because of this, we define u ja ( θ, γ ) asthe jump of log S a ( t ) at the point t j when it is assumed that the true crossing-time parameteris θ and the true initial dominance parameter is γ . More specifically, u ja ( θ, γ ) = log { S a ( t j ) } − log { S a ( t j − ) } , where S a ( t − ) denotes the left limit of S a ( t ) at time t . The term h ja ( θ, γ ) =1 − exp { u ja ( θ, γ ) } = 1 − S a ( t j ) /S a ( t j − ) can be interpreted as the discrete hazard at time point t j . Following Park and others (2012) and Johansen (1978), the log-likelihood function to bemaximized in this context islog L { u ( θ, γ ) , u ( θ, γ ) } = X a =0 m X j =1 h d ja log { − exp( u ja ( θ, γ )) } + ( R ja − d ja ) u ja ( θ, γ ) i , (2.5)where u ( θ, γ ) and u ( θ, γ ) are the m × u j ( θ, γ ) and u j ( θ, γ )respectively. In (2.5), R j = P ni =1 (1 − A i ) I ( Y i > t j ) denotes the number of individuals at riskin the control arm at time t j , and R j = P ni =1 A i I ( Y i > t j ) denotes the number of individualsat risk in the active treatment arm at time t j . In (2.5), d j = P ni =1 (1 − A i ) δ i I ( Y i = t j ) denotesthe number of events in the control arm at time t j while d j = P ni =1 A i δ i I ( Y i = t j ) denotes thenumber of events in the active treatment arm at time t j .In the absence of any crossing constraints on the survival functions, the only constraintson the vector u a ( θ, γ ) would be u ja ( θ, γ )
0, for j = 1 , . . . , m . Because the survival func-tions S ( t ) and S ( t ) associated with the vectors u ( θ, γ ) and u ( θ, γ ) have the form S a ( t ) =exp (cid:8) P mj =1 u ja ( θ, γ ) I ( t j t ) (cid:9) , the crossing constraints (2.3) and (2.4) can be expressed as acollection of linear inequality constraints. Specifically, we can represent both the crossing con-straints and the inequality constraints u ja ( θ, γ ) a k ( θ, γ )] T u ( θ, γ ) > k = 1 , . . . , m ,where a k ( θ, γ ) ∈ R m and u ( θ, γ ) = (cid:0) u ( θ, γ ) T , u ( θ, γ ) T (cid:1) T . For k m , the vectors a k ( θ, γ ) areconstructed to enforce the crossing constraints (2.3) or (2.4). If we define v ( θ ) = max { j : t j θ } if θ > t and v ( θ ) = 0 if θ < t , then a k ( θ,
1) for k m is given by a k ( θ,
1) = ((cid:0) k , m − k , − k , m − k (cid:1) , if k v ( θ ) (cid:0) − k , m − k , k , m − k (cid:1) , if v ( θ ) < k m where l denotes a vector of length l containing all ones, l denotes a vector of length l containingall zeros, and denotes a vector of “length 0” that should be ignored. Similarly, a k ( θ, −
1) is
N.C. Henderson & K. Nam & D. Feng given by a k ( θ, −
1) = ((cid:0) − k , m − k , k , m − k (cid:1) , if k v ( θ ) (cid:0) k , m − k , − k , m − k (cid:1) , if v ( θ ) < k m In order to enforce the constraints u ja ( θ, γ ) a k ( θ, γ ) is, for any value of ( θ, γ ), defined as a k ( θ, γ ) = ( k − m − , − , m − k ) for k = m + 1 , . . . , m .The maximum likelihood estimates ˆ u a ( θ, γ ) of u a ( θ, γ ) can be expressed as the solution tothe following optimization problem that has linear inequality constraintsmaximize log L (cid:8) u ( θ, γ ) , u ( θ, γ ) (cid:9) subject to A θ,γ (cid:20) u ( θ, γ ) T u ( θ, γ ) T (cid:21) > m , (2.6)where A θ,γ is the 3 m × m matrix whose k th row is a k ( θ, γ ) T . The above optimization prob-lem involves maximizing a concave function subject to linear inequality constraints, and hence,any local maximum is also guaranteed to be a global maximum (see e.g., Boyd and Vanden-berghe (2004)). In our implementation, we use sequential quadratic programming (Nocedal andWright (2006)) to compute the solution (cid:0) ˆ u ( θ, γ ) , ˆ u ( θ, γ ) (cid:1) of (2.6). Initialization of u ( θ, γ ), u ( θ, γ ) is done by minimizing, subject to the single-crossing constraint, the squared discrepancy P a =0 P mj =1 [ u ja ( θ, γ ) − log { ˆ S KMa ( t j ) / ˆ S KMa ( t j − ) } ] , where ˆ S KMa ( t ) is the Kaplan-Meier estimatesof the survival function in treatment arm a . One possible limitation of this computational strat-egy is that sequential quadratic programming can become very computationally demanding forlarge values of m . One remedy for this is to group the follow-up times Y , . . . , Y m into a collectionof small “bins” and set Y i to the midpoint of the bin to which it is assigned.2.3 Estimates of Crossing-Time Parameters θ and γ . The estimated vectors ˆ u ( θ, γ ) and ˆ u ( θ, γ ) will generate estimates of the two survival curves forfixed values of ( θ, γ ). To find the best values of the crossing parameters ( θ, γ ), we maximize theprofile log likelihood function ℓ P ( θ, γ ) associated with ˆ u ( θ, γ ) and ˆ u ( θ, γ ) ℓ P ( θ, γ ) = X a =0 m X j =1 h d ja log { − exp(ˆ u ja ( θ, γ )) } + ( R ja − d ja )ˆ u ja ( θ, γ ) i . (2.7)We refer to the values ˆ θ sc , ˆ γ sc which maximize ℓ P ( θ, γ ) as the single-crossing constrained estimatesof the crossing parameters θ and γ . Because the conditional estimates ˆ u a ( θ, γ ) do not change as θ varies over each of the intervals ( t j − , t j ) and can only change at each t j , the single-crossingconstrained estimates ˆ θ sc and ˆ γ sc can be found by solving the following discrete optimizationproblem (ˆ θ sc , ˆ γ sc ) = arg max θ ∈{ ,t ,...,t m − } ,γ ∈{− , } ℓ P ( θ, γ ). (2.8)The reason for only considering values of θ up to t m − is because both θ = 0 and θ = t m refer tosituations where one survival function dominates the other survival function at every time point t , . . . , t m . Thus, including θ = t m as a possible crossing time is superfluous as ( θ, γ ) = (0 ,
1) and( θ, γ ) = (0 , −
1) cover both scenarios where one survival curve dominates the other at each of theevent times t j .The vectors ˆ u a (ˆ θ sc , ˆ γ sc ) generate the following estimates of the survival functionsˆ S sca ( t ) = exp ( m X j =1 ˆ u ja (ˆ θ sc , ˆ γ sc ) I ( t j t ) ) . elayed Treatment Effects with Single-Crossing Constraints S sc ( t ) and ˆ S sc ( t ) as the single-crossing constrained estimates of the survival functions.Note that both ˆ S sc ( t ) and ˆ S sc ( t ) are flat for t > t m , and hence if ˆ S sc ( t ), ˆ S sc ( t ) satisfy the single-crossing constraint over [0 , t m ], they will satisfy it for all time points.2.4 Alternative Single-Crossing Constraints
The estimation strategy outlined in Sections 2.1-2.3 focuses on single-crossing constraints for thesurvival functions, but other related single-crossing constraints could potentially be incorporatedusing a similar approach. We briefly mention a few interesting possible extensions below. Whilewe explore non-smooth estimation of hazard functions with single-crossing constraints in ourapplication in Section 5, we do not explore the other mentioned extensions further as they liebeyond the scope of this paper.2.4.1
Non-smooth Estimation of Hazard Functions under Single-crossing Constraints
In manycases, it is more sensible to place single-crossing constraints on the hazards rather than on thesurvival curves. In this context, without imposing any smoothness conditions we would want thediscrete hazards h ja ( θ, γ ) = 1 − S a ( t j ) /S a ( t j − ) in one treatment arm to be larger (smaller) beforesome crossing time θ and remain smaller (larger) for j such that t j > θ . Because of the connection h ja ( θ, γ ) = 1 − exp { u ja ( θ, γ ) } between the u ja ( θ, γ ) and the discrete hazards h ja ( θ, γ ), we canexpress the single-crossing constraints on the discrete hazards as u j ( θ, γ ) > u j ( θ, γ ) for all j such that t j θu j ( θ, γ ) u j ( θ, γ ) for all j such that t j > θ, (2.9)if γ = 1 with both inequalities reversed whenever γ = −
1. As in the case of estimating survivalfunctions under single-crossing constraints, one would first, for fixed values of ( θ, γ ), find condi-tional maximum likelihood estimates of u ja ( θ, γ ) by maximizing the log-likelihood (2.5) subjectto constraints (2.9). After this, one would find estimates of the crossing parameters ( θ, γ ) bymaximizing the associated profile log-likelihood function (2.7).2.4.2 Smooth hazard functions with single-crossing constraints
To find smoothly-estimated haz-ard functions, one can consider hazard functions h a ( ·| θ, γ ) of the form h a ( t | θ, γ ) = 1 − exp n m X j =1 s ja ( t ) u ja ( θ, γ ) o = 1 − exp n s a ( t ) T u a ( θ, γ ) o , (2.10)for a choice of smoothing weights s a ( t ) = (cid:0) s a ( t ) , . . . , s m a a ( t ) (cid:1) T . A common choice of smooth-ing weights, for example, would be s ja ( t ) = b K { ( t − t j ) /b } for some symmetric kernel function K ( · ) and bandwidth b >
0. Under formulation (2.10), inequalities for the hazard functions attime points t j can be expressed as linear inequality constraints of the form s ( t j ) T u ( θ, γ ) > s ( t j ) T u ( θ, γ ), and hence, to compute smooth estimates of the hazard functions one could esti-mate the ˆ u ja ( θ, γ ) by maximizing the log-likelihood function (2.5) subject to the linear inequalityconstraints implied by the form of the hazard functions in (2.10). One would then find estimatesof the crossing-time parameters ( θ, γ ) by maximizing the associated profile log-likelihood function(2.7).An alternative to this approach would be to simply smooth estimated discrete hazards thathave been found using the approach outlined in Section 2.4.1. While this may work well in many N.C. Henderson & K. Nam & D. Feng situations, this approach would not guarantee that the smoothed hazard function estimates willsatisfy the single-crossing constraint.2.4.3
Covariate Adjustment
Suppose each individual in the study has an additional covariatewhich we denote with x i for the i th individual. In this context, we would let S a ( t | x i ) denote thesurvival function conditional on being assigned to treatment arm a and having covariate value x i and focus on crossing constraints for the “baseline” survival functions S ( t ) and S ( t ) where S a ( t ) = S a ( t | u ja ( θ, γ | x i ) = log { S a ( t j | x i ) /S a ( t j − | x i ) } and u ja ( θ, γ ) = u ja ( θ, γ | x i would be1 − exp { u ja ( θ, γ | x i ) } = (cid:2) − exp { u ja ( θ, γ ) } (cid:3) exp( β θ,γa x i ) , where β θ,γa is an arm-specific regression coefficient that can depend on the values of the crossingparameters. Note that this approach makes a proportional hazards assumption for the effect ofthe covariate x i within each treatment arm but does not make a proportional hazards assumptionfor the effect of the treatment arm assignment. Thus, this would still allow for the possibilityof having a crossing between the two baseline survival functions S ( t ) and S ( t ). If we assumedthat covariate-adjusted single-crossing constrained estimates of the baseline survival functionsmust have the same support as the overall single-crossing constrained estimates ˆ S sca ( t ) - namely,support on the observed event times t , . . . , t m - then to compute maximum likelihood estimatesof u ja ( θ, γ ) and ( β θ,γ , β θ,γ ) (for a fixed ( θ, γ )) in this context, we would maximize the followinglog-likelihood functionlog L cov { u ( θ, γ ) , u ( θ, γ ) , β θ,γ , β θ,γ } = X a =0 m X j =1 n X i =1 ˜ d jai h log { − exp( u ja ( θ, γ )) } + β θ,γa x i i + X a =0 m X j =1 n X i =1 ( ˜ R jai − ˜ d jai ) log n − (cid:2) − exp { u ja ( θ, γ ) } (cid:3) exp( β θ,γa x i ) o , where ˜ d jai = δ i I ( Y i = t j ) I ( A i = a ) and ˜ R jai = I ( Y i > t j ) I ( A i = a ). The assumption that thecovariate-adjusted single-crossing constrained estimates of S a ( t ) and the estimates ˆ S sca ( t ) havethe same support is similar to the assumptions made in, for example Owen (2001) and Zhou(2015), in the empirical likelihood analysis of the Cox proportional hazards model. Because therelationship between u ja ( θ, γ ) and the baseline survival functions in this context is the sameas the relationship between u ja ( θ, γ ) and the survival functions S a ( t ) in Sections 2.2-2.3, theconstraints on the u ja ( θ, γ ) required to impose a single-crossing constraint would be exactly thesame as those in (2.6).
3. Model Inference
Estimands of InterestMilestone Survival Probabilities.
Comparing differences in estimated survival probabilities atone or several pre-specified time points ˆ S sc ( t ∗ ) − ˆ S sc ( t ∗ ) , . . . , ˆ S sc ( t ∗ q ) − ˆ S sc ( t ∗ q ) can be a use-ful way of characterizing the treatment effect over time without relying on any assumptionsabout proportional hazards. Because ˆ θ sc provides an estimate of precisely where the sign changein S ( t ) − S ( t ) occurs, augmenting the survival probability differences at the milestones t ∗ , . . . , t ∗ q with the estimated crossing time ˆ θ sc can be helpful when interpreting the estimated differencesˆ S sc ( t ∗ j ) − ˆ S sc ( t ∗ j ). elayed Treatment Effects with Single-Crossing Constraints The Proportion Surviving up to Crossing.
The proportion surviving up to the crossing time intreatment arm a is represented by the parameter S a ( θ ). If both S and S are continuous, wewill have S ( θ ) = S ( θ ) and if either S or S are not continuous, these will be approximatelyequal as long as the true survival curves do not have large jumps. For this reason, we use S a ( θ ) todenote the proportion surviving up to crossing, and in practice, we estimate this parameter with { ˆ S sc (ˆ θ sc ) + ˆ S sc (ˆ θ sc ) } /
2. The quantity S a ( θ ) could be of particular interest if one is concernedabout a substantial fraction of patients experiencing early events that occur before the crossingtime. In these cases, reporting an estimate of S a ( θ ) provides a measure of the fraction of patientswho will survive long enough to reach the point at which the survival curve in the active treatmentarm begins to dominate to the control-arm survival curve. Restricted Mean Survival Time.
The restricted mean survival time (RMST) (see, e.g., Roystonand Parmar (2013)) is defined as the expected time under follow-up for an individual assumingyou only follow individuals up to some pre-specified time point τ . Specifically, the RMST forindividuals in treatment arm a is defined asRMST a ( τ ) = E { min( T i , τ ) | A i = a } .RMST a ( τ ) is equal to the area under the survival curve S a ( t ) between the time points t = 0 and t = τ , and the difference RMST ( τ ) − RMST ( τ ) provides an interpretable measure of treatmenteffect regardless of whether or not the proportional hazards assumption holds. While providing aninterpretable measure of treatment effect, the difference in RMST can mask important differencesin survival that occur at earlier time points. One way of addressing this is to also examinedifferences in RMST for different choices of τ with differences in RMST a ( θ ) perhaps being of keyinterest. Restricted Residual Mean Life.
If one is interested in differences in survival for those that arelonger survivors, the restricted residual mean life (RRML) function (see e.g., Cortese and others (2017)) is an appealing measure. The RRML function for treatment arm a is defined at time t asRRML a ( t, τ ) = E { min( T i , τ ) − t | A i = a, T i > t } = Z τt S a ( u ) S a ( t ) du .The quantity RRML a ( θ, τ ) represents the expected on-study survival conditional on the fact thatone has survived up to the crossing time θ . Crossing-Time Conditional Survival Curves.
In cases of delayed treatment where the two survivalcurves cross, it may be of interest to also plot survival probabilities conditional on surviving up tothe point of crossing. Such conditional probabilities give the probability of surviving past a pointof interest conditional on the fact that one has survived up to the crossing time. This conditionalsurvival curve for patients in treatment arm a is defined, for t > θ , as S a,cond ( t ) = P (cid:16) T i > t (cid:12)(cid:12)(cid:12) A i = a, T i > θ (cid:17) = S a ( t ) (cid:14) S a ( θ ). (3.11)The conditional survival curves may be estimated directly using ˆ S sca,cond ( t ) = ˆ S sca ( t ) / ˆ S sca (ˆ θ sc ). Itis worth mentioning that RRML a ( t, τ ) = R τθ S a,cond ( u ) du . Pre- and Post-crossing Average Hazard Ratios.
Comparing the average hazard ratios over thetime periods before and after the crossing can provide an interpretable measure of treatmentefficacy for longer survivors and can provide a good comparison for the relative improvement intreatment efficacy between earlier and later time points. Assuming arm-specific hazard functions0
N.C. Henderson & K. Nam & D. Feng h a ( t ) exist, we define, as in Kalbfleisch and Prentice (1981), the average hazard ratio using the“active treatment-to-total” hazard ratio which measures the average ratio between the activetreatment-arm hazard h ( t ) and the total hazard h ( t ) + h ( t ) across time. Specifically, for atruncation time τ and θ ∈ (0 , τ ), we define the pre- and and post-crossing average hazard ratiosas ¯ λ pre = 1 θ Z θ h ( t ) h ( t ) + h ( t ) dt and ¯ λ post = 1 τ − θ Z τθ h ( t ) h ( t ) + h ( t ) dt (3.12)respectively. One reason for using the ratio h ( t ) / { h ( t ) + h ( t ) } rather than h ( t ) /h ( t ) is toimprove estimation stability as potentially very small estimated value of h ( t ) could lead to highlyvariable estimates of ¯ λ pre and ¯ λ post . When assuming τ = t m , the parameters ¯ λ pre and ¯ λ post canbe estimated by the following quantitiesˆ¯ λ pre = 1ˆ θ sc m X j =1 ˆ h j (ˆ θ sc , ˆ γ sc )( t j − t j − ) I ( t j ˆ θ sc )ˆ h j (ˆ θ sc , ˆ γ sc ) + ˆ h j (ˆ θ sc , ˆ γ sc )ˆ¯ λ post = 1 t m − ˆ θ sc m X j =1 ˆ h j (ˆ θ sc , ˆ γ sc )( t j − t j − ) I ( t j > ˆ θ sc )ˆ h j (ˆ θ sc , ˆ γ sc ) + ˆ h j (ˆ θ sc , ˆ γ sc ) , where t = 0 and h ja ( θ, γ ) = 1 − exp { u ja ( θ, γ ) } is as defined in Section 2.2. Depending on thecontext, one could either use discrete hazard estimates ˆ h ja (ˆ θ sc , ˆ γ sc ) under the single-crossingconstraints on the survival functions or the single-crossing constraints on the hazard functiondescribed in Section 2.4.1. 3.2 Hypothesis Testing
While we can compute confidence intervals for certain parameters of interest that do not involvethe crossing parameters ( θ, γ ), it can be useful to perform inference with respect to both θ and another parameter of interest φ (or collection of parameters) that represents a measure oftreatment efficacy. For example, in traditional settings where it is assumed that proportionalhazards hold, φ would frequently be a hazard ratio, but in settings with delayed treatment effect,choosing φ to be an alternative estimand such as difference in RMST may be more appealing.Combining θ and φ in a joint hypothesis test can address concerns about having a scenario wherethe estimated value of φ indicates overall treatment effectiveness, but substantial time elapsesbefore the two survival curves clearly separate. Cases such as these may lead to concerns thatmost of the observed treatment benefit is mainly due to differences in long-term survivors.In the aforementioned context, one possible hypothesis of interest is that both the efficacyparameter φ is sufficiently large and the crossing time θ does not occur to late. This can beexpressed more formally as H : φ φ ∗ or θ > θ ∗ vs. H A : φ > φ ∗ and θ < θ ∗ , (3.13)where φ ∗ and θ ∗ are pre-specified values of ( θ, φ ) which are determined to be clinically meaningful.Alternatively, if it is difficult to specify a time point before which the crossing should occur,one could instead require that survival in the active treatment arm should be sufficiently largewhenever the crossing occurs. The hypothesis test of interest in this case would be H : φ φ ∗ or S ( θ ) p ∗ vs. H A : φ > φ ∗ and S ( θ ) > p ∗ . (3.14) elayed Treatment Effects with Single-Crossing Constraints φ and ˆ θ sc or ˆ φ and ˆ S sc (ˆ θ sc ) respectively. Another approach would be to test (3.13) or (3.14) by using abootstrap procedure to construct one-sided confidence intervals for either the parameter η =min { φ − φ ∗ , θ ∗ − θ } or the parameter η = min { φ − φ ∗ , S ( θ ) − p ∗ } . If the lower bound of theconfidence interval for η is greater than zero, one would reject H in (3.13). Likewise, a lowerbound for the confidence interval for η greater than zero would imply that one should reject H in (3.14).
4. Simulations
Estimation Performance with Piecewise Exponential Distributions
We considered six simulation scenarios where, in each scenario, it is assumed that survival followsa piecewise exponential distribution in both treatment arms. The arm-specific survival curves forthese six scenarios are depicted in Figure 2. The top-left graph in Figure 2 (Figure 2(a)) depictsScenario 1 where the survival curves never cross, and the survival curve for the active treatmentarm always dominates the control-arm survival curve. Figure 2 (b) depicts Scenario 2 where thereis a clear, unambiguous single crossing of the two survival curves at time point 5.0. Figure 2 (c)depicts Scenario 3 where the two survival curves have a single crossing near time point 2.0. WhileScenario 3 has a single, distinct crossing, when compared with Scenario 2 the two survival curvesin Scenario 3 do not have as much separation before the crossing occurs. Figure 2 (d) showsthe survival curves in Scenario 4 where there is a single crossing at time point 0.75, but in thisscenario, there is almost no separation between the curves before the crossing time. In Scenario5, there is a single crossing at time point 1.5 with relatively little separation before the crossingand a diminishing treatment benefit that occurs towards the end of the time interval considered.In Scenario 6, there are two crossing times, but the later crossing is more “distinct” than the firstin the sense that the separation between the two curves is larger immediately before and afterthe crossing point.For these simulations we set the total number of patients to n = 200, n = 400, and n = 800with the number of patients split evenly between the two treatment arms for each choice of n .For each of the six simulation scenarios and setting of n , we ran 200 simulation replications.The censoring distribution used in each of the six scenarios was a uniform distribution from 4to 8. While the percentage of survival outcomes which were observed event times varied acrosssimulation scenarios, the percentage was between 55% and 75% for each of the six settings, andthe percentage of observed events was typically 0 −
15% larger in the control arm than in theactive treatment arm. For each simulation setting, we evaluated the performance of the single-crossing constrained (SCC) procedure in estimating the following measures: difference in RMSTat time point 7, the differences in the survival function at the time points 2 and 4, the crossingtime θ , the proportion surviving up to crossing S a ( θ ), and the difference in RRML using the timepoints θ and 7.Table 1 shows the mean-squared error (MSE) for the SCC estimates across the six piecewise-exponential simulation scenarios and the three choices of sample size. For parameters that do notinvolve the crossing time θ , MSE was also computed for estimates based on the Kaplan-Meierestimates of the survival functions. As shown in this table, in scenarios with either no crossing ora distinct, single crossing (i.e., Scenarios 1-3), the SCC-based estimates had MSE performancewhich was consistently as good or better than the KM-based estimates for parameters for whichsuch a comparison could be made. For example, in the n = 400 settings, the reductions inMSE for the SCC-based estimator compared to the KM-based estimator of the RMST difference2 N.C. Henderson & K. Nam & D. Feng
Time S u r v i v a l P r obab ili t y (a) Scenario 1 Cross Time S u r v i v a l P r obab ili t y (b) Scenario 2 Cross Time S u r v i v a l P r obab ili t y (c) Scenario 3 Cross Time S u r v i v a l P r obab ili t y (d) Scenario 4 Cross Time S u r v i v a l P r obab ili t y (e) Scenario 5 Cross Cross Time S u r v i v a l P r obab ili t y (f) Scenario 6 Fig. 2. Arm-specific survival curves from the six piecewise exponential simulation scenarios considered inSection 4.
RMST (7) − RMST (7) were 2.3%, 3.5%, and 2.2% in Scenarios 1, 2, and 3 respectively, and inthe n = 800 settings the reductions in MSE for the RMST difference RMST (7) − RMST (7) were1.1%, 0.6%, and 0% in Scenarios 1, 2, and 3 respectively. MSE for the crossing-time estimatorˆ θ sc was lowest in both Scenarios 1 and 3 where there was either no crossing or an early, distinctcrossing. The relatively poorer result for ˆ θ sc in Scenario 2 is likely due to the fact that, inthis scenario, the true crossing time θ occurred much later in the study at a time where therewould typically be much fewer individuals remaining in this study. Despite this, the estimationperformance for the estimator of S a ( θ ) was quite good in Scenario 2 as both survival curvesare much more flat towards the end of the study period. Scenario 5 was the one setting whereestimation of the crossing-time parameter was notably poor. This was mainly due to the strongdiminished treatment effect present in Scenario 5 which often resulted in a crossing-time estimatecloser to the end of the considered time window rather than the much earlier true crossing timeof 1.5. Estimation performance of the SCC-based estimators were overall quite poor in Scenario elayed Treatment Effects with Single-Crossing Constraints ∆RMST(7) S (2) − S (2) S (4) − S (4) θ S a ( θ ) ∆RRML( θ, n Scenario SCC KM SCC KM SCC KM SCC SCC SCC200 1 0.1348 0.1362 0.0043 0.0044 0.0037 0.0036 0.6502 0.0221 0.15482 0.1254 0.1289 0.0041 0.0042 0.0047 0.0046 3.0721 0.0326 0.21283 0.1369 0.1368 0.0046 0.0044 0.0050 0.0049 1.3836 0.0182 0.16394 0.1540 0.1381 0.0054 0.0048 0.0055 0.0050 2.2295 0.0342 0.14725 0.1890 0.1265 0.0065 0.0045 0.0071 0.0053 7.4537 0.0728 0.31526 0.2047 0.1323 0.0068 0.0045 0.0064 0.0048 - - -400 1 0.0719 0.0736 0.0023 0.0024 0.0020 0.0020 0.3664 0.0092 0.07932 0.0669 0.0693 0.0019 0.0020 0.0025 0.0025 1.0848 0.0078 0.07283 0.0668 0.0683 0.0023 0.0023 0.0024 0.0024 0.6026 0.0088 0.07114 0.0730 0.0691 0.0024 0.0022 0.0025 0.0024 1.0826 0.0271 0.05955 0.1029 0.0636 0.0032 0.0022 0.0038 0.0025 7.4427 0.0699 0.19906 0.1184 0.0678 0.0035 0.0021 0.0036 0.0024 - - -800 1 0.0328 0.0332 0.0010 0.0010 0.0009 0.0009 0.0030 0.0008 0.03372 0.0317 0.0319 0.0009 0.0008 0.0011 0.0011 0.3767 0.0011 0.02493 0.0330 0.0330 0.0010 0.0010 0.0012 0.0012 0.2928 0.0043 0.03544 0.0360 0.0335 0.0012 0.0012 0.0012 0.0011 0.5016 0.0193 0.02785 0.0486 0.0299 0.0018 0.0012 0.0017 0.0011 5.3602 0.0459 0.11636 0.0610 0.0328 0.0022 0.0012 0.0017 0.0011 - - - Table 1. Mean-squared error (MSE) of single-crossing constrained (SCC) based estimators and Kaplan-Meier (KM) based estimators for several parameters from each of the six piecewise-exponential sim-ulation scenarios. MSE for the Kaplan-Meier based estimators are only shown when applicable. Theparameters ∆
RMST (7) and ∆
RRML ( θ,
7) are defined as ∆
RMST (7) = RMST (7) − RMST (7) and∆ RRML ( θ,
7) = RRML ( θ, − RRML ( θ,
5. Data Example
In this section, we examine reconstructed survival outcomes from a recently completed phase 3trial (Hellmann and others (2019)) examining the efficacy of a combination of immune checkpointinhibitors, nivolumab plus ipilimumab, for the treatment of non-small-cell lung cancer. In thistrial, patients were assigned to one of three treatments arms: a combination arm where nivolumabplus ipilimumab was administered, a monotherapy arm where nivolumab alone was administered,and a control arm where only chemotherapy was given. The primary endpoint in this studywas overall survival (OS) in the combination therapy arm versus the chemotherapy arm in thesubpopulation of patients whose tumors had an expression level of the programmed death ligand1 (PD-L1) that was at least 1%. Among the group of patients who had a PD-L1 expression of1% or more, 396 patients were assigned to the combination arm, and 397 patients were assignedto the chemotherapy only arm. While there was a notable delay in treatment effect in thisstudy, the analysis of this study reported in Hellmann and others (2019) concluded that thenivolumab plus ipilimumab treatment resulted in improved overall survival when compared withchemotherapy. In our analysis, we utilized survival outcomes that we reconstructed from thepublished Kaplan-Meier curves for OS in Hellmann and others (2019). Due to the resolution4
N.C. Henderson & K. Nam & D. Feng
Months S u r v i v a l P r obab ili t y Nivolumab + ipilimumabChemotherapy q ^ sc = 7.36 g ^ sc = 1 10 15 20 25 30 35 400.00.20.40.60.81.0 Months S u r v i v a l P r ob . Nivolumab + ipilimumabChemotherapy
Fig. 3. Single-crossing constrained estimates of survival in the reconstructed data from thenivolumab+iplimumab vs. chemotherapy trial. The left-hand panel shows the arm-specific survival curveestimates ˆ S sc ( t ) and ˆ S sc ( t ) along with estimates (ˆ θ sc , ˆ γ sc ) of the crossing-time parameters. The right-hand panel shows estimates ˆ S sc ,cond ( t ) and ˆ S sc ,cond ( t )of the crossing-time conditional survival curves; theconditional survival curves represent survival times conditional on surviving up to the crossing time. of these published images, our reconstructed survival outcomes are unlikely to be exactly thesame as those recorded in this study, but the reconstructed survival outcomes reproduce thepublished Kaplan-Meier curves quite closely. Using the reconstructed outcomes, median OS inthe nivolumab plus ipilimumab arm was 17.3 months while median OS in the chemotherapy armwas 15.0 months. While median OS suggests an overall benefit of the combination therapy, theKaplan-Meier estimates of OS indicate a delay in treatment effect as the estimated OS survivalcurve for the chemotherapy arm initially dominates the estimated OS survival curve for thecombination therapy arm, and a crossing appears to occur some time between 6 and 9 monthsbefore the two Kaplan-Meier estimates clearly separate at later time points.Figure 3 displays the single-crossing constrained estimates of the combination-arm and chemotherapy-arm survival curves for OS. As shown in this figure, the single-crossing constrained survival curveestimate for the chemotherapy arm shows an earlier superiority over the combination-arm sur-vival curve, while the combination-arm survival curve remains superior after the crossing occurs.The single-crossing constrained estimate ˆ θ sc of the crossing time was ˆ θ sc = 7.36 months, and thecorresponding estimate of the initial dominance parameter was ˆ γ sc = 1. The right-hand panelof Figure 3 shows the single-crossing constrained estimates of the conditional survival curves elayed Treatment Effects with Single-Crossing Constraints S a,cond ( t ) defined in (3.11). These curves represent estimates of survival probabilities conditionalon the fact that one has survived up to the crossing time. The graph of ˆ S sc ,cond ( t ) and ˆ S sc ,cond ( t )shows a clear superiority of the active treatment arm among those patients who will survive up toapproximately seven and a half months. Indeed, the probability for surviving more than 2 yearsconditional on surviving up to the crossing time is 0.54 in the combination arm and 0.47 in thechemotherapy arm, and the probability for surviving more than 3 years conditional on survivingup to the crossing is 0.44 in the combination arm and 0.29 in the chemotherapy arm.Table 2 displays single-crossing constrained estimates and their associated 95% confidenceintervals for other measures of treatment efficacy. To obtain these confidence intervals, we useda bootstrap with stratified resampling (Davison and Hinkley (1997)) where, in each bootstrapreplication, a subsample of the survival outcomes ( Y i , δ i ) was drawn with replacement from eachof the treatment arms. As shown in this table, our estimate of the proportion surviving up tocrossing parameter S a ( θ ) was 0.73 suggesting that approximately 73% of individuals in eithertreatment arm will survive up to the time point where the active treatment arm will begin tohave superior survival probabilities. The estimated difference in RMST truncated at 3 years was1.48 months. The estimated difference between the parameters RRML ( θ,
36) and RRML ( θ, θ S a ( θ ) 0.73 0.37 0.86RMST (36) − RMST (36) 1.48 -0.50 3.33RRML ( θ, − RRML ( θ,
36) 2.43 1.14 4.99 S (6) − S (6) -0.03 -0.09 0.02 S (12) − S (12) 0.04 -0.03 0.12 S (24) − S (24) 0.06 -0.01 0.13 S (36) − S (36) 0.11 0.05 0.18 S ,cond (12) − S ,cond (12) 0.06 0.00 0.15 S ,cond (24) − S ,cond (24) 0.08 0.00 0.16 S ,cond (36) − S ,cond (36) 0.15 0.08 0.24 Table 2. Single-crossing constrained estimates of different efficacy measures from the reconstructednivolumab+iplimumab vs. chemotherapy trial.
We also computed estimates of crossing-time parameters and the arm-specific hazard functionsunder a single-crossing constraint on the hazards rather than the survival curves. Here, we used theapproach described in Section 2.4.1 where a single-crossing constraint was placed on the discretehazards with the support of the discrete hazards being placed on the set of observed event times.The left-hand panel of Figure 4 shows the estimated discrete hazards for both treatment armswith the estimated hazard-crossing time of 2.4 months. This crossing-time estimate suggests that,while those in the combination arm initially have a larger hazard than those in the control arm,the advantage in hazard disappears roughly two and a half months before the hazards actuallycross. Using the crossing time of 2.4 months, estimates of the pre- and post-crossing averagehazard ratio parameters described in (3.12) were 0.77 and 0.32 respectively.While the hazard-based estimate of the crossing time can be useful, the discrete hazardsestimates are very non-smooth and hard to interpret. The right-hand panel of Figure 4 shows6
N.C. Henderson & K. Nam & D. Feng
Months ha z a r d Nivolumab + ipilimumabChemotherapy q ^ sc = 2.4 Months ha z a r d Nivolumab + ipilimumabChemotherapy
Fig. 4. Single-crossing constrained estimates of the hazard functions in the reconstructed data from thenivolumab+iplimumab vs. chemotherapy trial. The left-hand panel shows estimates of the arm-specificdiscrete hazards h ja ( θ, γ ) along with the estimate of 2.4 months for the crossing time for the discretehazards. The right-hand panel shows the smoothed hazard function for each treatment arm. hazard function estimates obtained by smoothing the discrete hazard estimates in the left-handpanel. To smooth the discrete hazards, we used the LOWESS smoother (Cleveland (1979)) withthe smoother span set to 2 /
3. We did not impose any additional single-crossing constraints whenperforming this smoothing, and for the time interval of 0 to 3 years, the single-crossing constraintfor the smoothed hazard functions was satisfied without requiring the use of additional constraintson the smoothed functions.
6. Conclusion
In this article, we have proposed nonparametric estimators of two survival curves when such curvesare constrained to cross at most once. The development of these single-crossing constrained esti-mators was primarily motivated by clinical trials involving recent cancer immunotherapies whereit is common to observe delays in treatment effect. While allowing for more than one crossingcould provide additional flexibility, our experience with immuno-oncology trials suggests thatmost successful therapies have at most one distinct crossing, and cases where one could arguethat multiple crossings are present in the underlying survival curves rarely provide clear evidenceof long-term benefit to patients. Though our approach can improve estimation performance incases where the underlying survival curves conform to a single-crossing constraint, one of the
EFERENCES S ( t ) and S ( t ). This could potentially improve power in caseswhere the active treatment shows a delayed treatment effect and where it is difficult to pre-specify the extent of the delay in treatment effect. One possible testing approach is to use aweighted Kaplan-Meier test statistic (Pepe and Fleming (1989)) with a weight function that isonly positive at time points after the estimated crossing time. This would closely resemble thetest statistic proposed by Logan and others (2008) who used a pre-specified rather than estimatedtime point to determine the support of their weight function. Another attractive alternative wouldbe to consider a weighted log-rank test with a piecewise constant weight function similar to theone proposed in Xu and others (2017) and, using the single-crossing constraint on the hazards,specify the jump of the weight function to occur at the estimated crossing time of the hazards.Though establishing the asymptotic null distribution of either the weighted Kaplan-Meier orweighted log-rank test statistic may be challenging, Monte Carlo permutation tests could beused to estimate the desired p-values. Even if the single-crossing constrained estimates are notused in constructing a test for comparing the arm-specific survival curves, another use of thesingle-crossing constrained estimates of the crossing parameters is in the design stage of a study.If relevant historical data are available, one could compute estimates of the crossing time of thehazards or survival functions and such estimates could be used to better inform parameter choicesused in sample size and power calculations. Supplemental Information An R package DelayedSurvFit implementing the methods described in this article and containingthe reconstructed dataset analyzed in Section 5 is publicly available at https://github.com/nchenderson/DelayedSurvFit . The R code used to conduct the simulation study described inSection 4 and the R code used for the data analysis shown in Section 5 are available at https://github.com/nchenderson/singlecrossingreproduce . ReferencesBoyd, Stephen and Vandenberghe, Lieven . (2004).
Convex optimization . Cambridge uni-versity press.
Chen, Tai-Tsang . (2013). Statistical issues and challenges in immuno-oncology.
Journal forimmunotherapy of cancer (1), 18. Chen, Tai-Tsang . (2015). Milestone survival: a potential intermediate endpoint for immunecheckpoint inhibitors.
Journal of the National Cancer Institute (9), djv156.8
REFERENCESCleveland, William S . (1979). Robust locally weighted regression and smoothing scatterplots.
Journal of the American statistical association (368), 829–836. Cortese, Giuliana, Holmboe, Stine A and Scheike, Thomas H . (2017). Regression modelsfor the restricted residual mean life for right-censored and left-truncated data.
Statistics inmedicine (11), 1803–1822. Davison, Anthony Christopher and Hinkley, David Victor . (1997).
Bootstrap methodsand their application , Number 1. Cambridge university press.
Demarqui, Fabio N and Mayrink, Vinicius D . (2019). A fully likelihood-based approach tomodel survival data with crossing survival curves. arXiv:1910.02406 . Harrington, David P and Fleming, Thomas R . (1982). A class of rank test procedures forcensored survival data.
Biometrika (3), 553–566. Hellmann, Matthew D, Paz-Ares, Luis, Bernabe Caro, Reyes, Zurawski, Bogdan,Kim, Sang-We, Carcereny Costa, Enric, Park, Keunchil, Alexandru, Aurelia,Lupinacci, Lorena, de la Mora Jimenez, Emmanuel and others . (2019). Nivolumab plusipilimumab in advanced non–small-cell lung cancer.
New England Journal of Medicine (21),2020–2031.
Johansen, Søren . (1978). The product limit estimator as maximum likelihood estimator.
Scandinavian Journal of Statistics (4), 195–199. Kalbfleisch, John D and Prentice, Ross L . (1981). Estimation of the average hazard ratio.
Biometrika (1), 105–112. Lawless, Jerald F . (2011).
Statistical models and methods for lifetime data , Volume 362. JohnWiley & Sons.
Lin, Ray S, Lin, Ji, Roychoudhury, Satrajit, Anderson, Keaven M, Hu, Tianle,Huang, Bo, Leon, Larry F, Liao, Jason JZ, Liu, Rong, Luo, Xiaodong and oth-ers . (2020). Alternative analysis methods for time to event endpoints under nonproportionalhazards: A comparative analysis.
Statistics in Biopharmaceutical Research (2), 187–198. Logan, Brent R, Klein, John P and Zhang, Mei-Jie . (2008). Comparing treatmentsin the presence of crossing survival curves: an application to bone marrow transplantation.
Biometrics (3), 733–740. Nocedal, Jorge and Wright, Stephen . (2006).
Numerical optimization . Springer Science& Business Media.
Owen, Art B . (2001).
Empirical likelihood . Chapman and Hall/CRC.
Pak, Kyongsun, Uno, Hajime, Kim, Dae Hyun, Tian, Lu, Kane, Robert C, Takeuchi,Masahiro, Fu, Haoda, Claggett, Brian and Wei, Lee-Jen . (2017). Interpretability ofcancer clinical trial results using restricted mean survival time as an alternative to the hazardratio.
JAMA oncology (12), 1692–1696. Park, Yongseok, Kalbfleisch, John D and Taylor, Jeremy MG . (2012). Constrainednonparametric maximum likelihood estimation of stochastically ordered survivor functions.
Canadian Journal of Statistics (1), 22–39. EFERENCES Pepe, Margaret Sullivan and Fleming, Thomas R . (1989). Weighted Kaplan-Meier statis-tics: a class of distance tests for censored survival data.
Biometrics , 497–507. Rahman, Rifaquat, Fell, Geoffrey, Ventz, Steffen, Arf´e, Andrea, Vanderbeek,Alyssa M, Trippa, Lorenzo and Alexander, Brian M . (2019). Deviation from theproportional hazards assumption in randomized phase 3 clinical trials in oncology: prevalence,associated factors, and implications.
Clinical Cancer Research (21), 6339–6345. Royston, Patrick and Parmar, Mahesh KB . (2013). Restricted mean survival time: analternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome.
BMC medical research methodology (1), 152. Schemper, Michael, Wakounig, Samo and Heinze, Georg . (2009). The estimation ofaverage hazard ratios by weighted Cox regression.
Statistics in medicine (19), 2473–2489. Xu, Zhenzhen, Zhen, Boguang, Park, Yongsoek and Zhu, Bin . (2017). Designing thera-peutic cancer vaccine trials with delayed treatment effect.
Statistics in medicine (4), 592–605. Yang, Song and Prentice, Ross . (2005). Semiparametric analysis of short-term and long-term hazard ratios with two-sample survival data.
Biometrika (1), 1–17. Zhao, Lihui, Claggett, Brian, Tian, Lu, Uno, Hajime, Pfeffer, Marc A, Solomon,Scott D, Trippa, Lorenzo and Wei, LJ . (2016). On the restricted mean survival timecurve in survival analysis.
Biometrics (1), 215–221. Zhou, Mai . (2015).