Nonparametric Analysis of Delayed Treatment Effects using Single-Crossing Constraints

Clinical trials involving novel immuno-oncology (IO) therapies frequently exhibit survival profiles which violate the proportional hazards assumption due to a delay in treatment effect, and in such settings, the survival curves in the two treatment arms may have a crossing before the two curves eventually separate. To flexibly model such scenarios, we describe a nonparametric approach for estimating the treatment arm-specific survival functions which constrains these two survival functions to cross at most once without making any additional assumptions about how the survival curves are related. A main advantage of our approach is that it provides an estimate of a crossing time if such a crossing exists, and moreover, our method generates interpretable measures of treatment benefit including crossing-conditional survival probabilities and crossing-conditional estimates of restricted residual mean life. We demonstrate the use and effectiveness of our approach with a large simulation study and an analysis of reconstructed outcomes from a recent combination-therapy trial.


Introduction
Recent advances in immuno-oncology (IO) therapies for the treatment of cancer have led to the development of a variety of treatments which show great potential for improving long-term patient outcomes.While very promising, patient response to such immunotherapies is often quite different when compared to more traditional cytotoxic agents such as chemotherapy.Indeed, it is well-recognized that IO drugs frequently exhibit a clear delay in treatment effect when they are compared with standard chemotherapies, and moreover, the nature of the delay in this treatment effect is often such that the estimated survival curves in the two treatment arms have a crossing at some time point after randomization.Due to this feature of immunotherapies, traditional comparisons 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 conventional choice for testing the superiority or non-inferiority of an active treatment over a control, and whenever the log-rank test passes a threshold for statistical significance, a hazard ratio is typically reported as one of the chief measures of treatment efficacy.However, in settings with survival curve crossings where the proportional hazards assumption is plainly violated, the interpretation of both the log-rank test and an estimated hazard ratio from a Cox proportional hazards model could be unclear.A variety of alternatives to the hazard ratio have been suggested in the context of evaluating the efficacy of immunotherapies or other context where the proportional hazards assumption 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 that are functionally "cured" (Chen (2013)).In the context of testing the equality of the arm-specific survival curves under possible delayed treatment effects, a number of weighted log-rank tests have been suggested including, for example, the piecewise weighted log-rank test (APPLE) proposed in Xu and others (2017) and weighted log-rank tests from the the Fleming-Harrington family which can place more weight on later time points (Harrington and Fleming (1982), Rahman and others (2019)).Combination tests which combine multiple weighted log-rank test statistics have also been proposed for analyzing IO trials (Lin and others (2020)).
While the measures of treatment effect listed above have a clear interpretation in the absence of proportional hazards and reporting such measures of treatment efficacy can certainly be valuable, imposing estimation constraints on how the survival curves may cross can provide additional information by which to evaluate a treatment that exhibits delayed benefit.To this end, we propose nonparametric estimates of the treatment-specific survival functions that allow for at most one crossing of the two survival functions without requiring that the crossing time be pre-specified.The approach we describe for estimating the survival curves under a single-crossing constraint can be thought of as a two-stage procedure.In the first stage, one finds conditional estimates of the survival function by maximizing a nonparametric log-likelihood that is conditional on the value of two crossing parameters.Then, in the second stage, one estimates these crossing parameters by maximizing a profile log-likelihood function.This produces estimates of the two survival curves which satisfy the single-crossing constraint, and it generates an estimate of the crossing time and an estimate of which survival curve is initially dominant.While flexible semiparametric approaches allowing for the crossing of treatment-specific survival functions have 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 related other than that they cross at most once, and moreover, our approach provides a likelihood-based estimate 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 there is a delay in treatment effect and improved long-term survival in the active treatment arm are often better modeled as having a single crossing after which point the survival curve in the active treatment arm consistently dominates the survival curve of the control arm.Such patterns of delayed 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 and provide greater interpretability when comparing the two arm-specific estimates of the survival functions.An important advantage of using single-crossing constraints is that it yields estimates and uncertainty intervals for the time that the crossing occurs.The estimate of the crossing-time parameter can be useful in both assessing when the active treatment begins to show superiority and as an interpretable measure that can be used as a component of additional measures of treatment efficacy that are designed for scenarios with delayed treatment effect which we outline in detail in Section 3.These efficacy measures include the proportion of patients surviving up to crossing, crossing-time conditional restricted residual mean life, crossing-time conditional survival probabilities, and pre/post-crossing average hazard ratios.
This paper is organized as follows.In Section 2, we describe a two-stage nonparametric estimation procedure which assumes that the treatment arm-specific survival curves cross at most once but otherwise makes no additional assumptions about the forms of the survival functions, and we briefly outline several potential extensions of this approach.Section 3 discusses a number of estimands of interest that can capture relevant concerns about the treatment impact in cases where survival curves may cross, and we describe how our method can be used to estimates these terms.Section 3 also describes hypothesis tests of interest under delayed treatment effect.Section 4 shows the results from a simulation study which evaluates the estimation performance of our method across six piecewise-exponential simulation settings.Section 5 shows an application of our method to a recent trial involving a novel combination immunotherapy in the treatment of non-small-cell lung cancer, and Section 6 concludes with a brief discussion.

Survival Curve Profiles and Notation
We assume that n patients have been enrolled in a randomized clinical trial consisting of two treatment 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 of directly observing T i , we observe the follow-up time Y i = min{T i , C i } and the event indicator 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 censoring is noninformative in the sense that the terms from the censoring distribution can be factored out of the likelihood function (Lawless (2011)), and we assume that C i and the treatment arm assignment A i are independent.
We let S a (t) = P (T i > t|A i = a) denote the survival function for those patients assigned to treatment group a ∈ {0, 1}.We consider an analysis where these survival functions are allowed to exhibit four distinct "survival profiles" according to whether or not a survival curve crossing occurs.Specifically, we allow for the possibility that the two survival functions S 0 (t) and S 1 (t) "cross", but we limit the number of crossings so that they can occur at most once.The crossing restriction implies that we will have one of four survival profiles, where each profile refers to a distinct crossing pattern of the survival functions.These four possible survival profiles are as follows: (1) a survival profile where treatment arm a = 1 completely dominates a = 0, namely S 1 (t) S 0 (t) for all t 0; (2) a survival profile where a = 0 dominates a = 1 before some crossing time θ but a = 1 dominates a = 0 afterwards, i.e., S 0 (t) S 1 (t) for 0 t θ and S 1 (t) S 0 (t) for t > θ; (3) a survival profile where a = 1 dominates a = 0 before some crossing time θ but a = 0 dominates a = 1 afterwards; and (4) a survival profile where a = 0 completely dominates 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 crossing parameters θ 0 and γ ∈ {−1, 1}.The parameter θ represents the crossing time of the survival functions, and the discrete parameter γ determines which treatment arm has the initially dominant survival function.Specifically, when γ = 1 the survival functions are assumed to exhibit the following behavior S 0 (t) S 1 (t) for 0 t θ and S 0 (t) S 1 (t) for t > θ, (2.1) and when γ = −1, the survival functions are assumed to obey the following inequalities S 0 (t) S 1 (t) for 0 t θ and S 0 (t) S 1 (t) for t > θ. (2.2) If there is a θ > 0 such that either constraints (2.1) or constraints (2.2) hold, we will say that θ is a non-trivial crossing time of the survival functions S 0 and S 1 .If there is no such θ > 0, we The crossing parameters θ and γ refer to the crossing-time parameter 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 γ = −1 indicates that the active treatment arm is dominant before the crossing.
must either have S 0 (t) S 1 (t) or S 1 (t) S 0 (t) for all t 0 in which case we set θ to the trivial crossing time θ = 0.When constructing estimates of the survival functions, we can only enforce the infinitedimensional constraints implied by (2.1) and (2.2) at a finite number of time points.In our implementation, we enforce the survival function constraints (2.1) and (2.2) at each of the observed event times.To this end, we let 0 < t 1 < t 2 < ... < t m denote the unique, ordered event times from both treatment arms, and hence, m represents the total number of events occurring in either of the treatment arms.For γ = 1 and a fixed value of θ, we enforce the following constraints S 0 (t j ) S 1 (t j ) for all j such that t j θ S 0 (t j ) S 1 (t j ) for all j such that t j > θ. (2.3) Similarly, for γ = −1 and a fixed value of θ, we enforce the following constraints S 0 (t j ) S 1 (t j ) for all j such that t j θ S 0 (t j ) S 1 (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 each choice of (θ, γ).In the next subsection, we describe our procedure for estimating the arm-specific survival functions when both θ and γ are assumed to be fixed.

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 constrained nonparametric maximum likelihood estimate of survival functions under a stochastic ordering constraint, the constrained maximum likelihood estimates of the survival functions will be discrete with potential jumps only at the event times t 1 , . . ., t m .Because of this, we define u ja (θ, γ) as the jump of log S a (t) at the point t j when it is assumed that the true crossing-time parameter is θ and the true initial dominance parameter is γ.More specifically, ) 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 be maximized in this context is where u 0 (θ, γ) and u 1 (θ, γ) are the m × 1 vectors having components u j0 (θ, γ) and u j1 (θ, γ) respectively.In (2.5), ) denotes the number of individuals at risk in the control arm at time t j , and ) denotes the number of individuals at risk in the active treatment arm at time t j .In (2.5), ) denotes the number of events in the control arm at time t j while d j1 = n i=1 A i δ i I(Y i = t j ) denotes the number of events in the active treatment arm at time t j .
In the absence of any crossing constraints on the survival functions, the only constraints on the vector u a (θ, γ) would be u ja (θ, γ) 0, for j = 1, . . ., m.Because the survival functions S 0 (t) and S 1 (t) associated with the vectors u 0 (θ, γ) and u 1 (θ, γ) have the form S a (t) = exp m j=1 u ja (θ, γ)I(t j t) , the crossing constraints (2.3) and (2.4) can be expressed as a collection of linear inequality constraints.Specifically, we can represent both the crossing constraints and the inequality constraints where 1 l denotes a vector of length l containing all ones, 0 l denotes a vector of length l containing all zeros, and 0 0 denotes a vector of "length 0" that should be ignored.Similarly, a k (θ, −1) is given by In order to enforce the constraints u ja (θ, γ) 0, a k (θ, γ) is, for any value of (θ, γ), defined as The maximum likelihood estimates ûa (θ, γ) of u a (θ, γ) can be expressed as the solution to the following optimization problem that has linear inequality constraints where A θ,γ is the 3m × 2m matrix whose k th row is a k (θ, γ) T .The above optimization problem 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 Vandenberghe (2004)).In our implementation, we use sequential quadratic programming (Nocedal and Wright (2006)) to compute the solution û0 (θ, γ), û1 (θ, γ) of (2.6).Initialization of u 0 (θ, γ), u 1 (θ, γ) is done by minimizing, subject to the single-crossing constraint, the squared discrepancy , where ŜKM a (t) is the Kaplan-Meier estimates of the survival function in treatment arm a.One possible limitation of this computational strategy is that sequential quadratic programming can become very computationally demanding for large values of m.One remedy for this is to group the follow-up times Y 1 , . . ., Y m into a collection of small "bins" and set Y i to the midpoint of the bin to which it is assigned.

Estimates of Crossing-Time Parameters θ and γ.
The estimated vectors û0 (θ, γ) and û1 (θ, γ) will generate estimates of the two survival curves for fixed values of (θ, γ).To find the best values of the crossing parameters (θ, γ), we maximize the profile log likelihood function ℓ P (θ, γ) associated with û0 (θ, γ) and û1 (θ, γ) (2.7) We refer to the values θsc , γsc which maximize ℓ P (θ, γ) as the single-crossing constrained estimates of the crossing parameters θ and γ.Because the conditional estimates ûa (θ, γ) do not change as θ varies over each of the intervals (t j−1 , t j ) and can only change at each t j , the single-crossing constrained estimates θsc and γsc can be found by solving the following discrete optimization problem ( θsc , γsc ) = arg max θ∈{0,t1,...,tm−1},γ∈{−1,1} ℓ P (θ, γ). (2.8) The reason for only considering values of θ up to t m−1 is because both θ = 0 and θ = t m refer to situations where one survival function dominates the other survival function at every time point t 1 , . . ., 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 the event times t j .
We refer to Ŝsc 0 (t) and Ŝsc 1 (t) as the single-crossing constrained estimates of the survival functions.Note that both Ŝsc 0 (t) and Ŝsc 1 (t) are flat for t t m , and hence if Ŝsc 0 (t), Ŝsc 1 (t) satisfy the singlecrossing constraint over [0, t m ], they will satisfy it for all time points.

Alternative Single-Crossing Constraints
The estimation strategy outlined in Sections 2.1-2.3 focuses on single-crossing constraints for the survival functions, but other related single-crossing constraints could potentially be incorporated using a similar approach.We briefly mention a few interesting possible extensions below.While we explore non-smooth estimation of hazard functions with single-crossing constraints in our application in Section 5, we do not explore the other mentioned extensions further as they lie beyond the scope of this paper.

Non-smooth Estimation of Hazard Functions under Single-crossing Constraints
In many cases, it is more sensible to place single-crossing constraints on the hazards rather than on the survival curves.In this context, without imposing any smoothness conditions we would want the discrete hazards h ja (θ, γ) = 1−S a (t j )/S a (t j −) in one treatment arm to be larger (smaller) before some 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 can express the single-crossing constraints on the discrete hazards as for all j such that t j θ u j0 (θ, γ) u j1 (θ, γ) for all j such that t j > θ, (2.9) if γ = 1 with both inequalities reversed whenever γ = −1.As in the case of estimating survival functions under single-crossing constraints, one would first, for fixed values of (θ, γ), find conditional maximum likelihood estimates of u ja (θ, γ) by maximizing the log-likelihood (2.5) subject to constraints (2.9).After this, one would find estimates of the crossing parameters (θ, γ) by maximizing the associated profile log-likelihood function (2.7).

Smooth hazard functions with single-crossing constraints
To find smoothly-estimated hazard functions, one can consider hazard functions h a (•|θ, γ) of the form for a choice of smoothing weights s a (t) = s 1a (t), . . ., s maa (t) T .A common choice of smoothing weights, for example, would be s ja (t) = 1 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 at time points t j can be expressed as linear inequality constraints of the form s 1 (t j ) T u 1 (θ, γ) s 0 (t j ) T u 0 (θ, γ), and hence, to compute smooth estimates of the hazard functions one could estimate the ûja (θ, γ) by maximizing the log-likelihood function (2.5) subject to the linear inequality constraints implied by the form of the hazard functions in (2.10).One would then find estimates of 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 that have been found using the approach outlined in Section 2.4.1.While this may work well in many situations, this approach would not guarantee that the smoothed hazard function estimates will satisfy the single-crossing constraint.
2.4.3Covariate Adjustment Suppose each individual in the study has an additional covariate which we denote with x i for the i th individual.In this context, we would let S a (t|x i ) denote the survival 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 0 (t) and S 1 (t) where S a (t) = S a (t|0).If we define u ja (θ, γ|x i ) = log{S a (t j |x i )/S a (t j − |x i )} and u ja (θ, γ) = u ja (θ, γ|0), a proportional hazards assumption with respect to the covariate x i would be where β θ,γ a is an arm-specific regression coefficient that can depend on the values of the crossing parameters.Note that this approach makes a proportional hazards assumption for the effect of the covariate x i within each treatment arm but does not make a proportional hazards assumption for the effect of the treatment arm assignment.Thus, this would still allow for the possibility of having a crossing between the two baseline survival functions S 0 (t) and S 1 (t).If we assumed that covariate-adjusted single-crossing constrained estimates of the baseline survival functions must have the same support as the overall single-crossing constrained estimates Ŝsc a (t) -namely, support on the observed event times t 1 , . . ., t m -then to compute maximum likelihood estimates of u ja (θ, γ) and (β θ,γ 0 , β θ,γ 1 ) (for a fixed (θ, γ)) in this context, we would maximize the following log-likelihood function where djai = The assumption that the covariate-adjusted single-crossing constrained estimates of S a (t) and the estimates Ŝsc a (t) have the 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 the relationship between u ja (θ, γ) and the baseline survival functions in this context is the same as the relationship between u ja (θ, γ) and the survival functions S a (t) in Sections 2.2-2.3, the constraints on the u ja (θ, γ) required to impose a single-crossing constraint would be exactly the same as those in (2.6).

Estimands of Interest
Milestone Survival Probabilities.Comparing differences in estimated survival probabilities at one or several pre-specified time points Ŝsc ) can be a useful way of characterizing the treatment effect over time without relying on any assumptions about proportional hazards.Because θsc provides an estimate of precisely where the sign change in S 1 (t)−S 0 (t) occurs, augmenting the survival probability differences at the milestones t * 1 , . . ., t * q with the estimated crossing time θsc can be helpful when interpreting the estimated differences Ŝsc 1 (t * j ) − Ŝsc 0 (t * j ).
The Proportion Surviving up to Crossing.The proportion surviving up to the crossing time in treatment arm a is represented by the parameter S a (θ).If both S 1 and S 0 are continuous, we will have S 1 (θ) = S 0 (θ) and if either S 1 or S 0 are not continuous, these will be approximately equal as long as the true survival curves do not have large jumps.For this reason, we use S a (θ) to denote the proportion surviving up to crossing, and in practice, we estimate this parameter with { Ŝsc 1 ( θsc ) + Ŝsc 0 ( θsc )}/2.The quantity S a (θ) could be of particular interest if one is concerned about a substantial fraction of patients experiencing early events that occur before the crossing time.In these cases, reporting an estimate of S a (θ) provides a measure of the fraction of patients who will survive long enough to reach the point at which the survival curve in the active treatment arm begins to dominate to the control-arm survival curve.
Restricted Mean Survival Time.The restricted mean survival time (RMST) (see, e.g., Royston and Parmar (2013)) is defined as the expected time under follow-up for an individual assuming you only follow individuals up to some pre-specified time point τ .Specifically, the RMST for individuals in treatment arm a is defined as 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 1 (τ ) − RMST 0 (τ ) provides an interpretable measure of treatment effect regardless of whether or not the proportional hazards assumption holds.While providing an interpretable measure of treatment effect, the difference in RMST can mask important differences in survival that occur at earlier time points.One way of addressing this is to also examine differences in RMST for different choices of τ with differences in RMST a (θ) perhaps being of key interest.
Restricted Residual Mean Life.If one is interested in differences in survival for those that are longer 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 as du.
The quantity RRML a (θ, τ ) represents the expected on-study survival conditional on the fact that one has survived up to the crossing time θ.
Crossing-Time Conditional Survival Curves.In cases of delayed treatment where the two survival curves cross, it may be of interest to also plot survival probabilities conditional on surviving up to the point of crossing.Such conditional probabilities give the probability of surviving past a point of interest conditional on the fact that one has survived up to the crossing time.This conditional survival curve for patients in treatment arm a is defined, for t > θ, as The conditional survival curves may be estimated directly using Ŝsc a,cond (t) = Ŝsc a (t)/ Ŝsc a ( θsc ).It is worth mentioning that RRML a (t, τ ) = τ θ S a,cond (u)du.Pre-and Post-crossing Average Hazard Ratios.Comparing the average hazard ratios over the time periods before and after the crossing can provide an interpretable measure of treatment efficacy for longer survivors and can provide a good comparison for the relative improvement in treatment efficacy between earlier and later time points.Assuming arm-specific hazard functions 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 active treatment-arm hazard h 1 (t) and the total hazard h 0 (t) + h 1 (t) across time.Specifically, for a truncation time τ and θ ∈ (0, τ ), we define the pre-and and post-crossing average hazard ratios as respectively.One reason for using the ratio h 1 (t)/{h 0 (t) + h 1 (t)} rather than h 1 (t)/h 0 (t) is to improve estimation stability as potentially very small estimated value of h 0 (t) could lead to highly variable estimates of λpre and λpost .When assuming τ = t m , the parameters λpre and λpost can be estimated by the following quantities λpre = 1 θsc where t 0 = 0 and h ja (θ, γ) = 1 − exp{u ja (θ, γ)} is as defined in Section 2.2.Depending on the context, one could either use discrete hazard estimates ĥja ( θsc , γsc ) under the single-crossing constraints on the survival functions or the single-crossing constraints on the hazard function described in Section 2.4.1.

Hypothesis Testing
While we can compute confidence intervals for certain parameters of interest that do not involve the 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 of treatment efficacy.For example, in traditional settings where it is assumed that proportional hazards 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 where the estimated value of φ indicates overall treatment effectiveness, but substantial time elapses before the two survival curves clearly separate.Cases such as these may lead to concerns that most 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 efficacy parameter φ is sufficiently large and the crossing time θ does not occur to late.This can be expressed more formally as (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 large whenever the crossing occurs.The hypothesis test of interest in this case would be One could test either (3.13) or (3.14) using a permutation test with the test statistics φ and θsc or φ and Ŝsc 1 ( θsc ) respectively.Another approach would be to test (3.13) or (3.14) by using a bootstrap procedure to construct one-sided confidence intervals for either the parameter η 1 = min{φ − φ * , θ * − θ} or the parameter η 2 = min{φ − φ * , S 1 (θ) − p * }.If the lower bound of the confidence interval for η 1 is greater than zero, one would reject H 0 in (3.13).Likewise, a lower bound for the confidence interval for η 2 greater than zero would imply that one should reject H 0 in (3.14).

Estimation Performance with Piecewise Exponential Distributions
We considered six simulation scenarios where, in each scenario, it is assumed that survival follows a piecewise exponential distribution in both treatment arms.The arm-specific survival curves for these six scenarios are depicted in Figure 2. The top-left graph in Figure 2 (Figure 2(a)) depicts Scenario 1 where the survival curves never cross, and the survival curve for the active treatment arm always dominates the control-arm survival curve.Figure 2 (b) depicts Scenario 2 where there is 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.While Scenario 3 has a single, distinct crossing, when compared with Scenario 2 the two survival curves in Scenario 3 do not have as much separation before the crossing occurs.Figure 2 (d) shows the survival curves in Scenario 4 where there is a single crossing at time point 0.75, but in this scenario, there is almost no separation between the curves before the crossing time.In Scenario 5, there is a single crossing at time point 1.5 with relatively little separation before the crossing and 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 first in the sense that the separation between the two curves is larger immediately before and after the crossing point.
For these simulations we set the total number of patients to n = 200, n = 400, and n = 800 with 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 4 to 8.While the percentage of survival outcomes which were observed event times varied across simulation scenarios, the percentage was between 55% and 75% for each of the six settings, and the percentage of observed events was typically 0 − 15% larger in the control arm than in the active treatment arm.For each simulation setting, we evaluated the performance of the singlecrossing constrained (SCC) procedure in estimating the following measures: difference in RMST at time point 7, the differences in the survival function at the time points 2 and 4, the crossing time θ, the proportion surviving up to crossing S a (θ), and the difference in RRML using the time points θ and 7.
Table 1 shows the mean-squared error (MSE) for the SCC estimates across the six piecewiseexponential simulation scenarios and the three choices of sample size.For parameters that do not involve the crossing time θ, MSE was also computed for estimates based on the Kaplan-Meier estimates of the survival functions.As shown in this table, in scenarios with either no crossing or a distinct, single crossing (i.e., Scenarios 1-3), the SCC-based estimates had MSE performance which was consistently as good or better than the KM-based estimates for parameters for which such a comparison could be made.For example, in the n = 400 settings, the reductions in MSE for the SCC-based estimator compared to the KM-based estimator of the RMST difference RMST 1 (7) − RMST 0 (7) were 2.3%, 3.5%, and 2.2% in Scenarios 1, 2, and 3 respectively, and in the n = 800 settings the reductions in MSE for the RMST difference RMST 1 (7)−RMST 0 (7) were 1.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, distinct crossing.The relatively poorer result for θsc in Scenario 2 is likely due to the fact that, in this scenario, the true crossing time θ occurred much later in the study at a time where there would typically be much fewer individuals remaining in this study.Despite this, the estimation performance for the estimator of S a (θ) was quite good in Scenario 2 as both survival curves are much more flat towards the end of the study period.Scenario 5 was the one setting where estimation of the crossing-time parameter was notably poor.This was mainly due to the strong diminished treatment effect present in Scenario 5 which often resulted in a crossing-time estimate closer to the end of the considered time window rather than the much earlier true crossing time of 1.5.Estimation performance of the SCC-based estimators were overall quite poor in Scenario 6, but this was a scenario where the assumption of a single crossing was plainly violated.

Data Example
In this section, we examine reconstructed survival outcomes from a recently completed phase 3 trial (Hellmann and others (2019)) examining the efficacy of a combination of immune checkpoint inhibitors, nivolumab plus ipilimumab, for the treatment of non-small-cell lung cancer.In this trial, patients were assigned to one of three treatments arms: a combination arm where nivolumab plus 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 study was overall survival (OS) in the combination therapy arm versus the chemotherapy arm in the subpopulation of patients whose tumors had an expression level of the programmed death ligand 1 (PD-L1) that was at least 1%.Among the group of patients who had a PD-L1 expression of 1% or more, 396 patients were assigned to the combination arm, and 397 patients were assigned to the chemotherapy only arm.While there was a notable delay in treatment effect in this study, the analysis of this study reported in Hellmann and others (2019) concluded that the nivolumab plus ipilimumab treatment resulted in improved overall survival when compared with chemotherapy.In our analysis, we utilized survival outcomes that we reconstructed from the published Kaplan-Meier curves for OS in Hellmann and others (2019).Due to the resolution of these published images, our reconstructed survival outcomes are unlikely to be exactly the same as those recorded in this study, but the reconstructed survival outcomes reproduce the published Kaplan-Meier curves quite closely.Using the reconstructed outcomes, median OS in the nivolumab plus ipilimumab arm was 17.3 months while median OS in the chemotherapy arm was 15.0 months.While median OS suggests an overall benefit of the combination therapy, the Kaplan-Meier estimates of OS indicate a delay in treatment effect as the estimated OS survival curve for the chemotherapy arm initially dominates the estimated OS survival curve for the combination therapy arm, and a crossing appears to occur some time between 6 and 9 months before the two Kaplan-Meier estimates clearly separate at later time points.
Figure 3 displays the single-crossing constrained estimates of the combination-arm and chemotherapyarm survival curves for OS.As shown in this figure, the single-crossing constrained survival curve estimate for the chemotherapy arm shows an earlier superiority over the combination-arm survival 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 the corresponding estimate of the initial dominance parameter was γsc = 1.The right-hand panel of Figure 3 shows the single-crossing constrained estimates of the conditional survival curves S a,cond (t) defined in (3.11).These curves represent estimates of survival probabilities conditional on the fact that one has survived up to the crossing time.The graph of Ŝsc 0,cond (t) and Ŝsc 1,cond (t) shows a clear superiority of the active treatment arm among those patients who will survive up to approximately seven and a half months.Indeed, the probability for surviving more than 2 years conditional on surviving up to the crossing time is 0.54 in the combination arm and 0.47 in the chemotherapy arm, and the probability for surviving more than 3 years conditional on surviving up 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% confidence intervals for other measures of treatment efficacy.To obtain these confidence intervals, we used a bootstrap with stratified resampling (Davison and Hinkley (1997)) where, in each bootstrap replication, a subsample of the survival outcomes (Y i , δ i ) was drawn with replacement from each of the treatment arms.As shown in this table, our estimate of the proportion surviving up to crossing parameter S a (θ) was 0.73 suggesting that approximately 73% of individuals in either treatment arm will survive up to the time point where the active treatment arm will begin to have superior survival probabilities.The estimated difference in RMST truncated at 3 years was 1.48 months.The estimated difference between the parameters RRML 1 (θ, 36) and RRML 0 (θ, 36) was 2.43 months which indicates that, conditional on surviving up to the crossing time, the expected gain in survival time was roughly two and a half months over the time period which begins at the crossing time and ends at 3 years.We also computed estimates of crossing-time parameters and the arm-specific hazard functions under a single-crossing constraint on the hazards rather than the survival curves.Here, we used the approach described in Section 2.4.1 where a single-crossing constraint was placed on the discrete hazards 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 arms with 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 actually cross.Using the crossing time of 2.4 months, estimates of the pre-and post-crossing average hazard ratio parameters described in (3.12) were 0.77 and 0.32 respectively.

Parameter
While the hazard-based estimate of the crossing time can be useful, the discrete hazards estimates are very non-smooth and hard to interpret.The right-hand panel of Figure 4 shows hazard function estimates obtained by smoothing the discrete hazard estimates in the left-hand panel.To smooth the discrete hazards, we used the LOWESS smoother (Cleveland (1979)) with the smoother span set to 2/3.We did not impose any additional single-crossing constraints when performing this smoothing, and for the time interval of 0 to 3 years, the single-crossing constraint for the smoothed hazard functions was satisfied without requiring the use of additional constraints on the smoothed functions.

Conclusion
In this article, we have proposed nonparametric estimators of two survival curves when such curves are constrained to cross at most once.The development of these single-crossing constrained estimators was primarily motivated by clinical trials involving recent cancer immunotherapies where it is common to observe delays in treatment effect.While allowing for more than one crossing could provide additional flexibility, our experience with immuno-oncology trials suggests that most successful therapies have at most one distinct crossing, and cases where one could argue that multiple crossings are present in the underlying survival curves rarely provide clear evidence of long-term benefit to patients.Though our approach can improve estimation performance in cases where the underlying survival curves conform to a single-crossing constraint, one of the

Fig. 1 .
Fig.1.Four possible survival profiles.The crossing parameters θ and γ refer to the crossing-time parameter 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 γ = −1 indicates that the active treatment arm is dominant before the crossing.

Fig. 2 .
Fig. 2. Arm-specific survival curves from the six piecewise exponential simulation scenarios considered in Section 4.

Fig. 3 .
Fig. 3. Single-crossing constrained estimates of survival in the reconstructed data from the nivolumab+iplimumab vs. chemotherapy trial.The left-hand panel shows the arm-specific survival curve estimates Ŝsc 1 (t) and Ŝsc 0 (t) along with estimates ( θsc, γsc) of the crossing-time parameters.The righthand panel shows estimates Ŝsc1,cond (t) and Ŝsc 0,cond (t)of the crossing-time conditional survival curves; the conditional survival curves represent survival times conditional on surviving up to the crossing time.

Fig. 4 .
Fig. 4. Single-crossing constrained estimates of the hazard functions in the reconstructed data from the nivolumab+iplimumab vs. chemotherapy trial.The left-hand panel shows estimates of the arm-specific discrete hazards hja(θ, γ) along with the estimate of 2.4 months for the crossing time for the discrete hazards.The right-hand panel shows the smoothed hazard function for each treatment arm.

Table 2 .
Single-crossing constrained estimates of different efficacy measures from the reconstructed nivolumab+iplimumab vs. chemotherapy trial.