A bounded hazard ratio Cox model for the effect of time to treatment on mortality

Abstract

In resource-limited or time-sensitive care settings, there is interest in assessing the impact of time to treatment (TTT) on mortality. Traditional Cox proportional hazards models, which specify the effect of TTT as an unrestricted term in the log hazard ratio, can produce counterintuitive results where the survival probability may not decrease monotonically with longer delays. Moreover, hazard ratios from such models quantify the effect of TTT conditional on surviving until treatment, rather than the effect of delayed treatment at baseline. We propose a class of bounded hazard ratio (BHR) Cox models that constrain the hazard ratio for TTT to attenuate towards the null with increasing treatment delay, such that hazard for death after treatment cannot exceed the hazard without treatment. Estimation can be performed using direct optimization of the partial log-likelihood or with an iterative linearized estimation procedure for large sample sizes. From BHR models, the estimated hazard ratio curve describes how treatment benefit diminishes with delay. Additionally, we propose a survival probability difference that provides an absolute measure for comparing survival under different treatment delays. We evaluate model performance in simulations and apply the method to examine treatment delays for colon cancer with data from the National Cancer Database.

Keywords

Time-to-event time-to-treatment survival semi-competing risks proportional hazards

1. Introduction

The effect of time to treatment (TTT) on subsequent survival outcomes has been examined in various clinical applications, such as cancer care,^1,2 emergency medicine,^3,4 and organ transplantation.⁵ In resource-limited settings, such as during a pandemic or within a constrained healthcare system, a lack of understanding of how TTT affects mortality can lead to inefficient allocation of resources that results in adverse outcomes.^6,7 Understanding and quantifying the effect of treatment delays on mortality can provide decision-makers with the necessary tools to make decisions about delaying and triaging care.

In this article, we are interested in assessing the effect of the timing of an intermediate event (e.g. treatment) on a subsequent terminal event (e.g. death) in the presence of censoring of both events. We focus on extrinsic intermediate events, such as treatment, whose timing is modifiable.⁸ Our goal is to quantify the effect of TTT on survival from baseline, when a decision to delay treatment may be made. Common analytical approaches include Cox models with delayed entry or with a time-dependent treatment indicator,⁹ sequential Cox models or landmarking,^10,11 and multi-state models for the joint distribution of TTT and time to death.^12,13 These and other methods have been compared in the literature.^14–16

In many clinical contexts, it is biologically plausible that treatment is beneficial when delivered promptly, but that its benefit attenuates as delay increases. For example, in cancer care, tumor growth during the waiting period may allow disease progression or spread, reducing the effectiveness of subsequent treatment.^17,18 Similarly, delays in acute or progressive conditions may lead to irreversible physiological deterioration prior to intervention.^19,20 These considerations motivate modeling frameworks in which treatment benefit is greatest at baseline and decreases with increasing delay, while remaining no worse than no treatment over the follow-up period of interest. However, standard survival models do not enforce such behavior in the implied hazard functions.

Constrained estimation has been explored in survival analysis to reflect knowledge about hazard shapes or the direction of covariate effects. Shape restrictions have been imposed on baseline hazards, and order-restricted or monotone covariate effects have been developed within proportional hazards models.^21–24 These approaches constrain regression coefficients or regression functions, thus restricting the ordering of hazard ratios across covariate levels. However, they do not guarantee that the implied hazard ratio satisfies properties such as boundedness relative to a reference hazard or convergence toward a specified limit as a covariate increases. In contrast, we consider a modeling structure that constrains the post-treatment hazard relative to the untreated hazard, ensuring that treatment benefit attenuates monotonically with delay and does not reverse direction. To our knowledge, such boundedness constraints on hazard ratios have not been previously formalized within Cox models.

Our contribution in this article is twofold. First, we show that a proportional hazards formulation that leaves the treatment-delay effect unconstrained can produce implausible behavior, including non-monotonic relationships between treatment delay and mortality risk. Consequently, we propose a class of bounded hazard ratio (BHR) Cox models that constrain the post-treatment hazard to remain below the untreated hazard while allowing treatment benefit to attenuate with increasing delay. Second, we demonstrate a measure that quantifies the impact of treatment delay on survival based on differences in marginal survival probabilities that account for the possibility of death prior to treatment. We demonstrate the application of the proposed class of models and quantification measures for describing the effect of TTT among colon cancer patients, a setting where traditional analytic approaches have produced conflicting results.

The remainder of the article is structured as follows. Section 2 describes the current methods using a proportional hazards effect for TTT and their theoretical limitations. In Section 3, we introduce the BHR Cox model, which imposes a plausible restriction on the behavior of the post-treatment hazard as TTT increases. Section 4 describes an alternative measure of the TTT effect to augment the interpretation of the hazard ratio from the BHR model. In Section 5, we conduct a simulation study to evaluate the behavior of the proposed models in settings of varying strength of TTT effect, sample sizes, and rates of competing risk cumulative incidence. In Section 6, we apply the proposed models to assess the effect of delaying treatment for colon cancer patients. We conclude with a discussion of the advantages, limitations, and possible extensions of the proposed approach in Section 7.

2. Proportional hazards effect of time to treatment

Let $\tilde{R}$ , $\tilde{T}$ , and $C$ be the true TTT (intermediate event), time to death (terminal event), and time to censoring, respectively, measured from a common baseline (e.g. diagnosis). The observed times are denoted by $R = min (\tilde{R}, \tilde{T}, C)$ and $Y = min (\tilde{T}, C)$ with corresponding event indicators $δ_{R} = I (R = \tilde{R})$ and $δ_{T} = I (Y = \tilde{T})$ , for treatment and death, respectively. Let $X$ be a vector of baseline covariates. The observed data on each of $n$ individuals is $D = {r, y, δ_{R}, δ_{T}, X}$ .

2.1. Time-dependent Cox proportional hazards model

A common approach for modeling TTT effects is a Cox model with time-dependent covariates for the effects of treatment and time-to-treatment, given by

h (t | r, X) = h_{0} (t) \exp {η I_{{t \geq r}} + γ I_{{t \geq r}} r + β^{T} X},

(1)

where

h_{0} (t)

is the baseline hazard and a binary time-dependent covariate

I_{{t \geq r}}

indicates whether the subject has received treatment by time

t

, and a second, continuous time-dependent covariate

I_{{t \geq r}} r

takes the value of the TTT after treatment and 0 before treatment. This model is fit to data from individuals who do and do not receive treatment during follow-up. From this model, the hazards of death before and after treatment are given by

h_{02} (t | X) = h_{0} (t) \exp {β^{T} X} and h_{12} (t | r, X) = h_{0} (t) \exp {η + γ r + β^{T} X},

(2)

respectively. This model assumes that the hazard for the transition from treatment at time

r

to death is proportional to the hazard of death without treatment with proportionality constant

e^{η}

. Thus,

e^{η}

represents the effect of immediate treatment on hazard of death (i.e. when

r = 0

). The effect of TTT,

γ

, represents a linear effect for increasing TTT on the log hazard ratio. The effect of treatment on the log-hazard is then represented by

η + γ r

This formulation uses a time-dependent treatment indicator to preserve the appropriate risk-set structure, avoiding immortal time bias when estimating associations between treatment timing and mortality.⁹ The specification in equation (1) represents the effect of TTT through a linear term in $r$ , but alternative parameterizations are possible. For example, the effect of treatment delay may be represented using categorical indicators, spline functions, or other flexible transformations of $r$ . Stratified time-dependent Cox models or multi-state models may also be used to allow the baseline hazards before and after treatment to differ.²⁵ While such approaches increase flexibility in representing the association between treatment timing and mortality, they do not impose structural restrictions on the resulting hazard functions.

2.2. Limitations of an unrestricted time-to-treatment effect

When modeling the effect of TTT on mortality, it is often clinically plausible that treatment is beneficial when delivered promptly and that this benefit diminishes as treatment is delayed. Under such assumptions, the hazard of death after treatment should not exceed the hazard without treatment, and treatment benefit should attenuate monotonically as delay increases. These considerations motivate the following structural properties of the post-treatment hazard:

Beneficial treatment. The hazard of death at time $t$ should be greater before versus after treatment for all treatment times, that is, $h_{02} (t | X) \geq h_{12} (t | r, X)$ .

Monotone TTT effect. The hazard of death at time $t$ should increase with treatment time, that is, $h (t | r_{1}, X) \leq h (t | r_{2}, X)$ for $r_{1} < r_{2}$ . For survival times $r_{1} \leq t < r_{2}$ , the hazard of death given treatment at $r_{2}$ is equal to the pre-treatment hazard, that is, $h (t | r_{2}, X) = h_{02} (t | X)$ , which gives the previous condition.

Decaying TTT effect. With increasing TTT, the hazard of death after treatment should approach (but not exceed) the hazard before treatment, that is, $lim_{r \to \infty} h_{12} (t | r, X) = h_{02} (t | X)$ .

These properties reflect a structural assumption about the behavior of the treatment effect rather than a specific parametric form.

The Cox models presented in the previous section specify the total treatment effect on the log hazard of death as $η + γ r$ . We refer to $η$ as the immediate treatment effect and $γ$ as the linear time-to-treatment effect. The conditions above imply that treatment is never harmful, and to satisfy these restrictions on the hazard requires that $γ \leq - η / r$ , or equivalently that the multiplicative total treatment effect on the hazard either reduces or does not change the hazard of death, that is, $e^{η + γ r} \leq 1$ . However, this is not a restriction imposed in the standard Cox model. Thus, the overall effect of receiving treatment will become harmful when $γ$ is sufficiently large or at later treatment times for any fixed $γ$ .

In Figure 1, we demonstrate these behaviors in the setting of proportional hazards with no baseline covariates, and a protective treatment effect, that is, $h_{02} (t) = c$ , $h_{12} (t | r) = h_{02} (t) e^{η + γ r}, η < 0$ . Figure 1a shows that the hazard rate for death after treatment may be increased instead of decreased, implying that treatment becomes harmful if the delay in treatment is long enough. In Figure 1b, these effects are demonstrated on the marginal survival function (to be defined explicitly in Subsection 4), where we see that for increasing $γ$ the relationship between treatment time $r$ and survival probability at certain times $t$ is not monotonic, and the survival curves for later treatment times cross below those for earlier treatment times. Increasing the flexibility of the functional form for $r$ , for example through categorical coding or spline-based representations, does not resolve this issue because the model remains unconstrained and may still imply post-treatment hazards that exceed the untreated hazard for some values of $r$ . This limitation reflects the trade-off between unconstrained models that may be more flexible, but do not enforce structural constraints and may be misspecified when the true effect is monotonic. In contrast, the constrained formulation provides a stable and interpretable representation when attenuation is expected.

Figure 1.

Effects of specifying the effect of time to treatment (TTT) as a linear term in a proportional hazards model on (a) hazard rates and (b) survival functions, illustrating how increases in the linear TTT effect parameter $(γ)$ can produce non-monotonic hazards and potential reversals of the treatment effect as treatment delay increases. TTT $r \in {0, 2, 4, 6, 8}$ is shown as separate lines within each panel. Plots were generated assuming an exponential baseline hazard $h_{02} (t) = 0.2$ , with a post-treatment hazard of $h_{12} (t | r) = h_{02} (t) \exp (η + γ r)$ , an immediate treatment effect of $η = \log (0.5)$ , and $γ \in {0, 0.2, 0.4}$ across columns.

3. The bounded hazard ratio cox model

To resolve these contradictions in modeling the effect of TTT on mortality risk, we seek an alternative class of models for the effect of TTT, such that conditions 1 $-$ 3 above are satisfied. Under these conditions, we assume that the greatest treatment benefit occurs with immediate treatment and that this benefit diminishes monotonically with increasing TTT. Thus, to restrict the behavior of the hazard post-treatment and at later treatment times relative to the pre-treatment hazard, we propose a proportional hazards model for the effect of immediate treatment, incorporating the TTT effect through a parametric function that reduces the effect of immediate treatment towards the null with increasing TTT.

3.1. Model specification

We propose a BHR Cox model, specified as:

h (t | r, X) = h_{0} (t) {[e^{η} (1 - w (r, α)) + w (r, α)]}^{I_{{t \geq r}}} \exp {β^{'} X}

(3)

where

w (r, α) : [0, \infty) \mapsto [0, 1]

is defined to be a weight function, non-decreasing in

r

, with

w (0, α) = 0

and

lim_{r \to \infty} w (r, α) = 1

. In this model,

e^{η}

represents the hazard ratio for immediate treatment (

r = 0

). The weight function

w (r, α)

controls how the treatment effect attenuates with delay, interpolating between the hazard with immediate treatment and the hazard without treatment. In practice, the weight function may be any cdf defined on

[0, \infty)

, dependent on a parameter vector

α

. Two examples of weight functions are a log-logistic cdf given by

w (r; α_{1}, α_{2}) = {1 + \exp [- \frac{\log (r) - α_{1}}{α_{2}}]}^{- 1} .

(4)

w (r; α_{1}, α_{2}) = 1 - \exp (- α_{1} r^{α_{2}})

, which resembles the cdf of a Weibull random variable.

Under this specification, the post-treatment hazard ratio relative to the pre-treatment hazard is $e^{η} (1 - w (r, α)) + w (r, α)$ , which is bounded between $e^{η}$ (when $r = 0$ ) and 1 (as $r \to \infty$ ). This directly addresses the three properties listed in the previous subsection. First, because $w (r, α) \in [0, 1]$ , the post-treatment hazard is always less than or equal to the pre-treatment hazard when $e^{η} \leq 1$ , satisfying Property 1 (beneficial treatment). Second, because $w (r, α)$ is non-decreasing in $r$ , the hazard ratio increases with treatment delay, satisfying Property 2 (monotone TTT effect). Third, since $lim_{r \to \infty} w (r, α) = 1$ , the post-treatment hazard approaches the pre-treatment hazard as treatment delay increases, satisfying Property 3 (decaying TTT effect).

Figure 2 illustrates possible behaviors of the attenuating treatment effect with increasing treatment delays and the corresponding impact on survival probability. The plot of the hazard ratio demonstrates how the treatment effect diminishes as TTT increases for various choices of the weight function. Since the hazard ratio directly determines the relative hazard over time, these shapes translate into the corresponding survival curves that show a monotonic association between longer TTT and decreasing survival, reflecting the monotone attenuation imposed by the weight function. As TTT increases, the survival probability with treatment returns to that without treatment. Thus, unlike a linear effect of TTT on the log hazard ratio, the BHR model provides natural limits on the effect of treatment delays, preventing the implausible scenario in which delaying treatment would result in worse survival than not treating at all.

Figure 2.

Different shapes of weight function and their corresponding hazard ratio (left) and survival (right) curves. The right panel shows the effect on survival of treatment administered at four possible times as modulated by the shape of the hazard ratio depicted in the left panel. The survival probability without treatment is shown in the right panels in black.

3.2. Estimation

We desire to obtain estimates of $η, α, β$ in model (3) in a semi-parametric fashion, that is, under the proportional hazards model. The BHR model is estimated through the Cox partial likelihood, allowing left-truncated (delayed entry) data to be handled using the standard Cox risk-set formulation, in which individuals enter the risk set at their truncation time and contribute only thereafter. Because $η, α$ enter the hazard function nonlinearly, the model may no longer be fit using standard Cox software. However, the proportional hazards model may still be used to obtain estimates of both $η$ and the weight parameters $α$ , as the nonlinear form of the weight does not involve $t$ but only the time-dependent covariates $I_{{t \geq r}}$ and $r$ .

We adopt two approaches for estimation. First, we may directly optimize the partial likelihood function corresponding to this nonlinear-in-parameters hazard ratio using a general-purpose optimizer, such as Nelder–Mead. This approach is effective for moderate sample sizes, such as those considered in our simulation study ( $N \leq 10,000$ ), but may become prohibitively slow as the number of observations or covariates increases further.

As an alternative, we propose an iterative linearized estimation procedure based on linearizing the log hazard ratio and leveraging the efficient estimation available in coxph() from the survival package.²⁶ This procedure provides a fast and stable approach to estimating the partial likelihood maximizer by iteratively updating a locally linear surrogate model using standard Cox software while targeting the same partial likelihood as the original nonlinear specification. Estimation via this iterative method proceeds as follows.

Let $θ = (η, α^{⊤})^{⊤}$ denote the parameter vector governing the treatment effect and weight function, and let $g (θ) = \log [e^{η} (1 - w (r, α)) + w (r, α)]$ denote the log hazard ratio at treatment time $r$ . The estimation algorithm proceeds by repeatedly linearizing $g (θ)$ with respect to $θ$ and fitting the resulting surrogate model using standard Cox regression software, analogous to the iteratively reweighted least squares procedure used for generalized linear models. At iteration $k$ , a Taylor expansion gives

g (θ) \approx g ({\hat{θ}}^{(k)}) + \nabla g ({\hat{θ}}^{(k)})^{T} (θ - {\hat{θ}}^{(k)}) .

The gradient components are

\frac{\partial g}{\partial η} = \frac{e^{η} (1 - w (r, α))}{e^{η} (1 - w (r, α)) + w (r, α)}, \frac{\partial g}{\partial α} = \frac{(1 - e^{η}) \partial w (r, α) / \partial α}{e^{η} (1 - w (r, α)) + w (r, α)} .

This linear predictor can be fit as a Cox model with an offset and auxiliary covariates based on the gradient, and coefficients are updated until convergence. Baseline covariate effects

β

enter linearly and are re-estimated at each iteration via the Cox partial likelihood, and therefore do not influence the nonlinear component

g (θ)

. Starting values may be obtained from a pilot fit to a small subset of the data using the explicit partial likelihood for our nonlinear model. Standard errors may be obtained either directly from the final coxph object or from numerical approximations to the Hessian of the true log-likelihood function evaluated at the converged parameter estimates. The R code for estimation is available at https://github.com/ksuresh17/TTT-BHRCoxModel.

To test for the presence of a TTT effect in the form of the proposed BHR weight functions, a likelihood ratio test (LRT) may be used. However, since under the null hypothesis of no TTT effect the parameters of the weight function are either on the boundary of the parameter space²⁷ or not identified, we recommend a simulation-based approach to obtain the distribution of the LRT statistic under the null.²⁸ Details for this procedure are given in Supplementary Material B.3. For each simulated dataset, both the null model and the alternative BHR model (with either Weibull or log-logistic weight functions) are fit using the same estimation procedure as in the main analysis. The likelihood ratio statistic is then computed as $2 {\log L_{BHR} - \log L_{null}}$ for each replicate. The empirical distribution of these simulated LRT statistics may then be used as an estimate of the null distribution. The $p$ -values are calculated as the proportion of the null simulations with an LRT statistic more extreme than what was observed.

A subset of the parameter space may produce nearly identical weight functions over the observed treatment time range, leading to weak identifiability of the individual parameters in $α$ while preserving accurate estimation of the hazard ratio curve. The log-logistic and Weibull weight functions allow flexible attenuation patterns in the treatment effect as demonstrated in Figure 2, although alternative weight functions could also be specified if a particular clinical setting motivates a different attenuation pattern. To guide the selection of the weight function, one can compare model fits using standard likelihood-based criteria such as Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC).

4. Assessing the time-to-treatment effect

4.1. The bounded hazard ratio curve

By maintaining the proportional hazards framework, the BHR Cox models provide estimates of a hazard ratio curve representing the instantaneous risk of death at any time $t$ for a particular treatment time $r$ compared to remaining untreated up to time $t$ . The slope of the tangent line to the log hazard ratio curve at a given treatment time $r$ reflects how quickly the hazard is increasing with additional delay at that point. This quantity plays a similar conceptual role to the time-to-treatment coefficient $γ$ in a traditional Cox model, which assumes a constant rate of increase in log hazard. In contrast, the BHR model allows this slope to vary smoothly with $r$ , providing a more flexible representation of how treatment benefit attenuates over time.

The hazard ratio curve allows one to identify periods where the treatment effect is constant or minimal and treatment delays do not affect risk (flat portions of the curve), periods where there is a significant increase in risk for even small treatment delays (steep increases in the curve), and the threshold time after which receiving treatment provides a similar benefit to not receiving treatment at all (point at which the difference between the curve and its asymptote at 1 becomes negligible). Thus, these hazard ratio curves provide insight into optimizing the timing of interventions.

Figure 2 gives examples of three possible hazard ratio curve shapes from the same log-logistic family of weight functions. A concave curve may represent an emergency medicine setting such as sepsis where early administration of antibiotics is critical and each hour of delay increases mortality risk significantly; however, later in the clinical course the intervention timing is less crucial since the outcome is already determined. A convex curve may represent chronic diseases, such as type 2 diabetes or cardiovascular disease, where the timing of lifestyle interventions may be more lenient closer to diagnosis but becomes more important later in life as more comorbidities accumulate. A sigmoidal curve is a possible representation of the effect of treatment delay in certain cancers, such as breast or colon, where following early detection there are periods of slow followed by rapid tumor growth and a leveling off when the cancer reaches an advanced stage after which initiating treatment does not improve survival.

4.2. The survival probability difference

While the hazard ratio curve characterizes the instantaneous effect of treatment delay on mortality, it is inherently conditional on surviving until treatment. In many clinical settings, however, decisions about treatment timing are made at baseline (e.g. diagnosis), and delaying treatment also carries the risk of death prior to treatment. Thus, it is useful to consider a survival metric that quantifies the absolute survival probability associated with delaying treatment from baseline. We therefore consider the marginal survival probability under treatment at time $r$ , defined as

P (T > t | r, X) = I_{{t < r}} e^{- \int_{0}^{t} h_{02} (u | X) d u} + I_{{t \geq r}} e^{- \int_{0}^{r} h_{02} (u | X) d u} e^{- \int_{r}^{t} h_{12} (u | r, X) d u},

(5)

where

h_{02} (t | X)

and

h_{12} (t | r, X)

represent the hazard of death given covariates before and after treatment, respectively. Under the BHR Cox model in equation (3), this probability becomes

P (T > t | r, X) = I_{{t < r}} e^{- H_{0} (t) \exp {β^{T} X}} + I_{{t \geq r}} e^{- H_{0} (r) \exp {β^{T} X}} e^{- [H_{0} (t) - H_{0} (r)] \exp {β^{T} X} [e^{η} (1 - w (r, α)) + w (r, α)]},

where

H_{0} (t) = \int_{0}^{t} h_{0} (u) d u

is the cumulative baseline hazard. The first term represents survival prior to treatment and does not depend on the treatment time. The second term represents survival after treatment and incorporates the treatment effect. Although this probability conditions on treatment occurring at time

r

, it does not condition on surviving to receive treatment and therefore reflects the competing risk of death while waiting to receive treatment.

Notice that even when there is no time-to-treatment effect (i.e. $w (r, α) = 0$ ), the marginal survival probability will depend on the treatment time $r$ whenever there is an immediate treatment benefit ( $η \neq 0$ ) and a non-zero hazard of death prior to treatment. In this setting, earlier treatment still produces higher survival probabilities because delaying treatment increases the time during which individuals remain exposed to the pre-treatment mortality hazard. This illustrates an important distinction between the hazard ratio and marginal survival metrics, where a constant hazard ratio that does not vary with treatment delay may suggest no effect of TTT on mortality, while the SPD can still indicate improved survival with earlier treatment because deaths can occur while waiting for treatment. This situation may arise in acute clinical settings such as trauma or sepsis, where mortality risk is substantial before treatment is delivered.

We use this marginal survival probability to define a survival probability difference (SPD) comparing treatment at two times $r_{1} < r_{2}$ for an individual with covariate vector $X$ as

S P D_{r_{1}, r_{2}} (t | X) = P (T > t | r_{2}, X) - P (T > t | r_{1}, X),

(6)

where

P (T > t | r_{2}, X)

and

P (T > t | r_{1}, X)

are defined in equation (5). This quantity represents the change in survival probability at time

t

when treatment occurs at time

r_{2}

instead of

r_{1}

. Under the expectation that treatment is beneficial and

r_{2} > r_{1}

, it will generally hold that

P (T > t | r_{2}, X) < P (T > t | r_{1}, X)

, so that

S P D_{r_{1}, r_{2}} (t | X) < 0

. The SPD therefore quantifies the loss in survival probability associated with delaying treatment from

r_{1}

r_{2}

, and captures the cumulative impact of treatment timing on survival. In contrast to the hazard ratio, which describes the instantaneous relative risk conditional on survival, the SPD reflects differences in survival probabilities that accumulate over time. For ease of interpretation we will generally set

r_{1} = 0

, so that equation (6) represents the decrease in survival probability associated with treating at time

r_{2}

instead of immediately. Computation of SPD requires models that specify both the pre- and post-treatment hazards.

The magnitude of the SPD depends on both the pre-treatment hazard and how treatment effectiveness changes with delay. In settings where mortality risk is concentrated early after diagnosis or injury, even small treatment delays may produce large SPD values because delaying treatment prolongs exposure to the high pre-treatment hazard. Conversely, when early mortality risk is low and treatment effectiveness decays slowly with delay, delaying treatment may have little impact on survival probabilities. Alternative cumulative measures based on the marginal survival distribution, such as restricted mean survival time (RMST) differences, could also be considered.

5. Simulation study

5.1. Simulation design

We conduct a simulation study to assess the behavior of the proposed BHR Cox model under three different types of hazard ratio curves (concave, convex, sigmoidal) and three mortality cumulative incidence scenarios (20%, 40%, 75% cumulative incidence of the competing risk of death).

Individual data is simulated from the BHR Cox model in equation (3) assuming a Weibull baseline hazard for TTT and time to death, and a log-logistic cdf weight function, given in equation (4). Further details of data generation are given in Supplementary Material A.1. The true parameter values are given in Supplementary Tables A1 and A2, and correspond to the hazard ratio curves shown in Figure 2. Censoring times are simulated from an exponential distribution with rate parameters chosen to achieve 20 $-$ 30% censoring.

Due to the applicability of this method in large observational data sets, we compare the performance of the methods in sample sizes of $N = 2500$ , $5000$ , and $10 000$ . We additionally explore the estimation of the proposed BHR method in smaller samples of $500, 1000, 1500$ , and $2000$ individuals for all the hazard ratio curve shapes in the setting of 40% mortality cumulative incidence.

The slow attenuation of the hazard ratio curve in the convex setting was chosen to assess the performance of the model in settings with a weak time-to-treatment effect. We additionally consider a convex curve setting with a stronger time-to-treatment effect and faster attenuation to the null during the follow-up period. To demonstrate model performance when adjusting for covariates, we also consider a scenario (40% cumulative incidence of death, sample size 5000) that incorporates one binary and one continuous covariate associated with the hazard of death. We also assess performance under misspecification by generating data using a piecewise hazard ratio curve with an initial flat segment followed by a linear increase that does not follow the parametric BHR weight functions.

The data-generating model, denoted “BHR (log-logistic)”, as well as a misspecified BHR Cox model with a Weibull cdf weight function, denoted “BHR (Weibull)”, are fit to each simulated data set, with parameters estimated by optimizing the log partial likelihood function. We also estimate these models using the iterative linearization approach, denoted “BHR* (log-logistic)” and “BHR* (Weibull).” Since we have argued in Section 2 that the Cox model with a linear TTT effect on the log hazard ratio is conceptually inappropriate in this setting, we do not use that model as a data-generating mechanism; however, we do fit it as a comparator in the simulation study, denoted “Cox (linear),” as well as a Cox model with only the treatment effect but not a time-to-treatment predictor, denoted “Cox (null).” We additionally fit a flexible Cox model allowing the treatment effect to vary over time by interacting treatment with a natural cubic spline function of time (df = 4), denoted “Cox (spline).”

We simulate 500 data sets and assess estimation performance for the immediate treatment effect $η$ using bias, confidence interval coverage probability, empirical standard error (ESE), and asymptotic standard error (ASE). Estimation of the weight function parameters is evaluated using the root mean squared difference (RMSD) between the areas under the true and estimated hazard ratio curves, as their individual values may not be precisely recovered. RMSD summarizes integrated differences across the hazard ratio curve rather than at a particular delay time point. Bias and mean squared error (MSE) are also computed for the predicted survival probability difference (SPD) relative to immediate treatment at five treatment times selected to correspond to time points at which patients remain at risk and treatment occurs.

5.2. Results

Figure 3 shows estimated hazard ratio curves for the three hazard ratio curve shapes and varying sample sizes under the 40% mortality cumulative incidence scenario. Supplementary Figures A1–A3 display additional results for the 20%, 40%, and 75% scenarios, including interquartile ranges. Estimates for the immediate treatment effect and hazard ratio curves appear in Supplementary Tables A3–A5 and A6, respectively.

Figure 3.

Median of estimated hazard ratio curves with 40% mortality cumulative incidence for varying curve shapes and sample size (N). Results are presented for the direct optimization bounded hazard ratio models (BHR) with log-logistic and Weibull weight functions, and Cox models with no (null) and linear time to treatment effect. The data-generating model uses the log-logistic weight function (black solid line).

Bias in the estimates of $η$ is lower with decreasing cumulative incidence of mortality and larger sample sizes, and is highest under the concave hazard ratio curve (Supplementary Tables A3–A5). For concave hazard ratio curves, the treatment benefit attenuates rapidly after treatment initiation, so relatively little information is available about the immediate treatment effect beyond the earliest treatment times, which increases the difficulty of estimating $η$ . Root mean squared differences between the true and estimated areas under the hazard ratio curves also decrease with lower mortality incidence and increasing sample size, and are comparable across the three curve types (Supplementary Table A6). When including baseline covariates that are associated with the hazard of death, the parameter estimates from the BHR models remain unbiased and the BHR curve is well recovered (Supplementary Tables A7–A8).

Greater estimation error for the convex hazard ratio curve is likely due to its weak TTT effect, as the hazard ratio attenuates slowly to 1 in this setting (Figure 2). The misspecified Cox models also show low bias in $η$ and reasonable hazard curve estimates at early treatment times, due to the near-linear and weak TTT effect in this setting. In the setting of a steep convex hazard ratio curve, estimation is improved and is similar to that in the sigmoidal hazard ratio setting (Supplementary Material A.6). The Cox model with a spline function of TTT produced immediate treatment effect estimates that were similar to those from the BHR models across settings, although bias decreased more slowly with increasing sample size and empirical standard errors were larger. Since the Cox models do not enforce boundedness, both the linear and spline-based Cox models produce hazard ratio estimates that exceed 1 at later treatment times in the sigmoidal setting (Supplementary Figure A1), where the true hazard ratio attenuates toward the pre-treatment hazard. This implies a higher hazard of death with treatment than without treatment at longer delays, which is inconsistent with the assumed treatment benefit.

The BHR models estimated using the iterative linearized estimation procedure (BHR*) showed similar results for the bias of the estimated immediate treatment effect and hazard ratio curves as the models estimated with direct optimization (BHR) (Supplementary Tables A3–A5, A6). They did, however, suffer from some numerical issues for smaller sample sizes, especially in the setting of high mortality cumulative incidence. The median computation time for estimation of the BHR* models was lower than for BHR models across all hazard ratio curve shapes ( $N = 10,000$ : sigmoidal 1.1s vs. 2.3s; concave 1.6s vs. 113s; convex 1.5s vs. 119s). This reflects the computational advantage of the linearization approach, especially in settings where the likelihood surface is relatively flat.

Figure 4 shows bias in estimated survival probability differences (SPD) across hazard ratio curve types and mortality settings, using a sample size of 10,000. Supplementary Figures A4–A6 display corresponding results for other sample sizes, including interquartile ranges. We evaluate the bias in predicted SPD values when comparing delayed versus immediate treatment at five specified time points. The bias is low for the BHR models even for the smaller sample sizes and greater mortality cumulative incidence. In the concave setting, the SPD is underestimated by the misspecified BHR model and both the correctly specified and misspecified weight functions have greater variability in their predictions. In the sigmoidal and concave settings with stronger TTT effects, the Cox models overestimated the SPD, with the largest bias observed for the Cox models with null and linear TTT effects and attenuated bias for the model with spline effects. Reflecting the weak TTT effect in the convex setting, the bias of SPD is low for the misspecified Cox models.

Figure 4.

Median of the bias of the survival probability difference versus immediate treatment at survival time $t = 8$ with a sample size of 10,000 for varying hazard ratio curves and mortality cumulative incidence (MCI) settings. Results are presented for Cox models with no time to treatment effect (null) and a linear time to treatment effect (linear), and bounded hazard ratio models estimated using direct optimization (BHR) with a log-logistic and Weibull weight function. The data-generating model is the weighted Cox model with the log-logistic weight function (BHR (log-logistic)).

In Supplementary Material A.7, we demonstrate that for smaller sample sizes between 500 and 2000, the bias of the immediate treatment effect increases, with reduced convergence rates observed in some settings, consistent with weaker identification at small sample sizes (Supplementary Table A11). Estimation of the hazard ratio curve shape remains reasonably accurate at smaller sample sizes, with root mean squared differences (RMSD) generally below 0.05 for $N \geq 1500$ and improving with increasing $N$ (Supplementary Figure A9, Supplementary Figure A10).

To assess robustness to model misspecification, we simulated data under a piecewise hazard ratio curve that does not follow the parametric weight function used in the BHR model (Supplementary Material A.8). In this setting, the BHR models continued to provide accurate estimates of the immediate treatment effect and were able to recover the hazard ratio and SPD curves well. The Cox models with null and linear TTT effects produced biased estimates of the immediate treatment effect, while the spline-based model violated the expected monotonicity of the hazard ratio, exhibiting non-monotonic behavior in regions where the true hazard ratio was constant.

Overall, the BHR models, whether estimated via direct or iterative methods, exhibit low bias for the hazard ratio curves and survival probability differences (SPD) across a range of hazard curve shapes and sample sizes.

6. Colon cancer application

To examine the performance of the proposed methods in practice, we applied them to colon cancer data from the National Cancer Database (NCDB). Exploring the effect of TTT in colon cancer is of particular importance since treatment delays often occur during critical periods of high tumor growth rates¹ and have been associated with increased risk of metastasis and limiting treatment options.² Conflicting results have been reported for the association between TTT and survival in colon cancer, with several studies reporting negative or U-shaped associations between TTT and mortality.²⁹ While some of these inconsistencies may be attributed to differences in study population or adjustment variables, previous results are also often based on analytic methods that do not appropriately account for delayed entry or use time-dependent Cox models with an unrestricted hazard ratio.^6,18,30–33

We provide an analysis of this data that estimates the TTT effect within the proposed BHR framework. We compare the BHR model results with both a null model (no TTT effect) and a model with a linear TTT effect. Constraining the behavior of the hazard ratio to satisfy a conceptually sound model may result in more reasonable findings than the traditional Cox modeling approaches.

6.1. Data description and modeling approaches

The NCDB is a joint project of the Commission on Cancer of the American College of Surgeons and the American Cancer Society. We analyzed NCDB data on individuals diagnosed with colon cancer between 2010 and 2019, stratified by stage at diagnosis. Exclusion criteria and adjustment variables followed prior work⁶ and are detailed in Supplementary Material B.

Two standard Cox models were fit, one with only treatment included, “Cox (null)” and one with treatment plus its interaction with TTT, “Cox (linear).” Two BHR Cox models were fit, one with a Weibull cdf weight function, “BHR (Weibull),” and one with a log-logistic cdf weight function, “BHR (log-logistic).” These models differ in the structural assumptions imposed on how treatment benefit evolves with delay. Because of these differences, likelihood-based comparisons are not the primary basis for interpretation. Instead, we focus on the resulting hazard ratio and survival probability patterns. Expanded likelihood-based comparisons are provided in the supplement. Models were adjusted for potential demographic and clinical confounders, described in Supplementary Material B.

Additionally, as a diagnostic measure to assess calibration of the estimated survival probability differences, we compared the fit of these models to a more flexible Cox model with time-dependent strata for treatment time. In this model, there is no explicit term for TTT; instead the model allows a different baseline hazard for mortality in each TTT stratum. The first stratum was defined to be pre-treatment; then, patients were assigned to one of 10 strata based on a grid with intervals spaced 0.05 years apart based on their observed treatment time.

6.2. Estimated treatment and TTT effects

After exclusions, data was available on approximately $55 000$ patients diagnosed with stage I disease, $78 500$ with stage II, and $76 000$ with stage III. In this data, TTT is short relative to survival times, placing us in a setting with very low competing risk of death. Median TTT was 30 days for stage I and 26 days for stages II and III. Ninety percent of patients were treated within 77 days, 60 days, and 58 days for stages I, II, and III, respectively. Stage-specific Kaplan–Meier curves (Supplementary Material, Figure B1) are provided in the supplement as descriptive summaries. Expanded model comparisons, including $Δ$ log-likelihood and $Δ$ AIC across the Cox and BHR models, are also provided in Supplementary Material (Table B1).

Figure 5.

Hazard ratios for treatment as a function of time to treatment, as estimated from the colon cancer NCDB data. Shaded areas represent pointwise 95% confidence intervals for the hazard ratio. Note that the $y$ -axis is on the log scale.

Estimated hazard ratio curves from the four models for the effect of treatment on mortality are shown in Figure 5. These represent the ratio of the instantaneous risk of death for a patient treated at the time indicated on the $x$ -axis relative to remaining untreated, adjusted for important clinical and demographic factors. The value of the curve at TTT $= 0$ represents the effect of immediate treatment. Note that the null model shows a flat horizontal line, as this model’s hazard ratio is constant with respect to TTT. According to the likelihood ratio test with simulated null distribution, all BHR models fit the data significantly better than models with no TTT effect (all $p \leq 0.01$ ; Table 1).

Table 1.

Likelihood ratio test results comparing proposed BHR fits with a null model, that is, no TTT effect.

	Log-logistic		Weibull
Stage	LR	$p$ -Value	LR	$p$ -Value
Stage I	297.5	<0.001	292.2	<0.001
Stage II	418.1	<0.001	418.5	<0.001
Stage III	7.0	0.010	8.2	0.006

p-values obtained by simulating the null distribution of the test statistic.

There are notable differences in the estimated hazard ratio curves by stage at diagnosis. For stage I disease, the hazard ratio is less than 1 for all models and TTT $\leq 3$ months. This is consistent with what would be expected under a large beneficial treatment effect that is only weakly dependent on TTT. However, for patients diagnosed with stage II or III disease, the picture is less clear. Both the null and linear TTT-effect models estimate hazard ratios greater than 1, implying a reversal of the treatment effect. The BHR models, by contrast, show estimated hazard ratios of less than 1, with an estimated hazard ratio for immediate treatment of approximately 0.6 for stage II and 0.8 for stage III, respectively.

The BHR models suggest that the treatment effect then gradually returns to the null for treatment by 2–3 months post-diagnosis. Although these models produce different hazard ratio patterns, all are adjusted for the same covariates and remain subject to potential unobserved or time-varying confounding. Thus, the differences between models reflect the structural assumptions imposed on the treatment effect. The BHR model does not remove confounding bias, but instead estimates treatment effects under the assumption that treatment benefit attenuates smoothly toward the null as treatment delay increases.

6.3. Goodness of fit

Figure 6 shows the SPD distributions as median and upper and lower quartiles, at treatment times of 1, 2, 3, and 4 months. All SPD values were calculated with reference to immediate treatment, that is, $r_{1} = 0$ in equation (6), so values less than zero imply there is a survival benefit associated with earlier treatment. We observe that this is the case for all models other than the null (for stages II and III) and stratified models (at treatment time of 1 month). The BHR models closely track the stratified Cox model in stage I and show broadly similar trends in stages II and III, providing a constrained summary of the observed treatment-timing pattern. Because the Cox linear and BHR models are not nested and impose different structural assumptions on the hazard ratio, we compare their implied survival probability differences and their agreement with the stratified Cox calibration model rather than relying solely on likelihood-based criteria. The linear Cox model tends to overestimate the reduction in survival from delaying treatment, with SPD values that are more negative than those of the stratified model, particularly at longer delays. The BHR results align with those from the stratified Cox model while providing a parsimonious summary under the imposed monotone attenuation structure.

Figure 6.

SPD functions as estimated from the colon cancer NCDB data. These represent the difference in survival probability between hypothetical patients treated at the time indicated by the $x$ -axis value vs. immediate treatment; patients in this comparison are identical with respect to all other covariates. Points represent the median value across the sample while error bars represent the 25th and 75th percentiles of these quantities. A horizontal line at 0 is provided for reference.

7. Discussion

In resource-limited settings that require time-sensitive care, quantitative methods that characterize the effect of treatment delays are needed to guide prioritization and timing of care. In this article, we identify a mathematical inconsistency in modeling TTT as a predictor with an unconstrained effect on the hazard in a Cox model, which can imply that receiving treatment increases mortality risk. We introduce a BHR model that avoids this contradiction by incorporating the effect of immediate treatment and a parametric weight function that attenuates this effect with increasing delay.

Traditional Cox models that assume a log-linear treatment effect are not constrained to satisfy the structural assumptions considered here for the relationship between treatment delay and survival, as seen in our colon cancer application where they implied a harmful effect of immediate treatment. More flexible approaches, such as stratified or piecewise Cox models that estimate separate hazards across treatment-time categories, allow greater flexibility in the shape of the treatment-timing effect, but may produce non-monotonic patterns or be less stable in sparsely observed regions. The BHR model provides a compromise: it enforces smooth monotonic attenuation of treatment benefit while allowing the data to determine the rate of decay. This balance improves both stability and interpretability of the estimated hazard ratio curve across a range of admissible attenuation patterns, including near-linear, plateauing, and flat patterns.

The effect of treatment timing is often explored in large observational cohorts or registries with heterogeneous treatment timing, as in our colon cancer application. We propose an iterative linearized estimation procedure that leverages standard Cox model software for large samples where direct maximization of the partial likelihood is computationally prohibitive. The procedure builds on the gradient-based optimizer in coxph(), providing a fast and scalable alternative to direct optimization. For moderate sample sizes, direct maximization with Nelder–Mead was sufficiently robust and accurate. Although the method uses local linearization, it converges to the exact maximum partial likelihood solution and can be extended to other models with nonlinear parameter structures.

Simulation results showed that larger sample sizes are needed in settings with high rates of mortality cumulative incidence, an underlying concave hazard ratio curve, or modest TTT effects (i.e. slow attenuation of the hazard ratio curve to the null) where there is weaker identification of the weight function parameters from the data. In some settings, an approximately linear TTT effect from a traditional Cox model may be appropriate for the range of relevant treatment delay times; however, the BHR model remains preferable, as it can approximate linear effects while accommodating a broader range of plausible hazard ratio shapes.

Although our proposed BHR weight functions are chosen to represent realistic relationships between TTT and mortality risk, it is possible that they do not adequately capture the true shape of the TTT–mortality risk curve in some data sets. One can consider using more flexible weight functions, for example, the generalized or three-parameter gamma³⁴ cdf. These specifications, however, lack closed-form expressions for the derivatives of the weight function required by the proposed linearization estimation approach, and thus may have limited implementation in larger data sets.

Causal inference frameworks, such as marginal structural models and the clone-censor-weight approach, have been developed to estimate the effects of time-dependent treatments in the presence of time-dependent confounding.^35,36 Our approach does not address causal identification, but instead proposes a modeling framework for representing how treatment benefit evolves with delay while ensuring clinically plausible hazard behavior. In principle, the BHR formulation could be incorporated within causal inference frameworks, where the causal method addresses confounding while the BHR structure models how treatment benefit evolves over time.

There are other settings where a BHR Cox model may be appropriate. For example, in examining the effect of treatment discontinuation, one could use the proposed weight function to balance between the hazard of death pre- and post-discontinuation with the condition that the hazard on treatment never exceeds the hazard post-discontinuation. This can be implemented within the proposed framework with some reparameterization to account for a decreasing rather than increasing weight function, and can take advantage of the model flexibility to threshold at a non-null value if indicated by the data.

The BHR Cox model was proposed within the proportional hazards framework. Alternative survival models, such as an additive hazard model or a multi-state model that assumes non-proportional hazards for the transitions to death pre- and post-treatment also require imposing constraints to ensure the monotonicity of mortality risk with increasing TTT. Without the application of the weight function to a proportionality constant, this is not a direct extension of the proposed approach. Future research will explore how similar restrictions might be applied to investigate a TTT effect using other survival models.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802261458075 - Supplemental material for A bounded hazard ratio Cox model for the effect of time to treatment on mortality

Supplemental material, sj-pdf-1-smm-10.1177_09622802261458075 for A bounded hazard ratio Cox model for the effect of time to treatment on mortality by Krithika Suresh, Wenyu Zhang and John D Rice in Statistical Methods in Medical Research

Footnotes

Acknowledgments

The data used in the study are derived from a de-identified NCDB file. The American College of Surgeons and the Commission on Cancer have not verified and are not responsible for the analytic or statistical methodology employed, or the conclusions drawn from these data by the investigator.

ORCID iDs

Krithika Suresh

Wenyu Zhang

John D Rice

Ethical approval and informed consent statements

There are no human participants in this article, and informed consent was not required.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Research reported in this publication was supported by The Sunway Trust Scholarship and Internship funds and the National Cancer Institutes of Health under Award Number P30CA046592. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Declaration of conflicting interest

The authors have no conflicting interests to declare with respect to the research, authorship, and/or publication of this article.

Data availability statement

The data used in this study are available from the National Cancer Data Base (NCDB). Access to this data requires approval and can be requested directly from the NCDB.

Supplemental material

Supplemental material for this article is available online.

References

Amri

Bordeianou

Sylla

, et al.

Treatment delay in surgically-treated colon cancer: does it affect outcomes?

Ann Surg Oncol 2014; 21: 3909–3916.

Kaltenmeier

Shen

Medich

, et al. Time to surgery and colon cancer survival in the united states. Ann Surg 2021; 274: 1025–1031.

Band

Gaieski

Hylton

, et al. Arriving by emergency medical services improves time to treatment endpoints for patients with severe sepsis or septic shock. Acad Emerg Med 2011; 18: 934–940.

Suresh

Dixon

Patel

, et al. The epidemiology and outcomes of prolonged trauma care (epic) study: methodology of a prospective multicenter observational study in the western cape of south africa. Scand J Trauma Resusc Emerg Med 2022; 30: 55.

Meier-Kriesche

Port

Ojo

, et al. Effect of waiting time on renal transplant outcome. Kidney Int 2000; 58: 1311–1317.

Cone

Marchese

Paciotti

, et al. Assessment of time-to-treatment initiation and survival in a cohort of patients with common cancers. JAMA network open 2020; 3: e2030072.

Hartman

Sun

Devasia

, et al. Integrated survival estimates for cancer treatment delay among adults with cancer during the covid-19 pandemic. JAMA Oncol 2020; 6: 1881–1889.

Epstein

Muñoz

. A bivariate parametric model for survival and intermediate event times. Stat Med 1996; 15: 1171–1185.

Kalbfleisch

Prentice

. The statistical analysis of failure time data. Hoboken, NJ: J. Wiley, 2002.

10.

Andersen

Gill

. Cox’s regression model for counting processes: a large sample study. Ann Stat 1982; 10: 1100–1120.

11.

Van Houwelingen

Putter

. Dynamic prediction in clinical survival analysis. Boca Raton, FL: CRC Press, 2012.

12.

Meira-Machado

de Uña-Álvarez

Cadarso-Suárez

, et al. Multi-state models for the analysis of time-to-event data. Stat Methods Med Res 2009; 18: 195–222.

13.

Valeri

Proust-Lima

Fan

, et al. A multistate approach for the study of interventions on an intermediate time-to-event in health disparities research. Stat Methods Med Res 2023; 32: 1445–1460.

14.

Agarwal

Moshier

, et al. Immortal time bias in observational studies of time-to-event outcomes: assessing effects of postmastectomy radiation therapy using the national cancer database. Cancer Control 2018; 25: 1073274818789355.

15.

Jones

Fowler

. Immortal time bias in observational studies of time-to-event outcomes. J Crit Care 2016; 36: 195–199.

16.

Karim

Gustafson

Petkau

, et al. Comparison of statistical approaches for dealing with immortal time bias in drug effectiveness studies. Am J Epidemiol 2016; 184: 325–335.

17.

Hanna

King

Thibodeau

, et al. Mortality due to cancer treatment delay: systematic review and meta-analysis. bmj 2020; 371: m4087.

18.

Khorana

Tullio

Elson

, et al. Time to initial cancer treatment in the united states and association with survival over time: an observational study. PLoS ONE 2019; 14: e0213209.

19.

De Luca

Suryapranata

Ottervanger

, et al. Time delay to treatment and mortality in primary angioplasty for acute myocardial infarction. Circulation 2004; 109: 1223–1225.

20.

Kumar

Roberts

Wood

, et al. Duration of hypotension before initiation of effective antimicrobial therapy is the critical determinant of survival in human septic shock. Crit Care Med 2006; 34: 1589–1596.

21.

Chung

. Nonparametric estimation of proportional hazards with monotone baseline hazard and covariate effect. Stat Biosci 2024; 16: 787–800.

22.

Chung

Ivanova

Hudgens

, et al. Partial likelihood estimation of isotonic proportional hazards models. Biometrika 2018; 105: 133–148.

23.

Dunson

Herring

. Bayesian inferences in the cox model for order-restricted hypotheses. Biometrics 2003; 59: 916–923.

24.

Gelfand

Mallick

. Bayesian analysis of proportional hazards models built from monotone functions. Biometrics 1995; 51: 843–852.

25.

Le-Rademacher

Therneau

. The utility of multistate models: a flexible framework for time-to-event data. Curr Epidemiol Rep 2022; 9: 183–189.

26.

Therneau

. A Package for Survival Analysis in R. 2022. https://CRAN.R-project.org/package=survival. R package version 3.4-0.

27.

Self

Liang

. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. J Am Stat Assoc 1987; 82: 605–610.

28.

Rice

Johnson

Juarez-Colunga

, et al. Application of gap time analysis with flexible hazards to pulmonary exacerbations in the epic observational study. Biomet J 2022; 64: 1075–1089.

29.

Franssen

Strous

Bongers

, et al. The association between treatment interval and survival in patients with colon or rectal cancer: a systematic review. World J Surg 2021; 45: 2924–2937.

30.

Ang

Lee

Chia

, et al. Time-to-treatment initiation and its effect on all-cause mortality: insights from the surveillance, epidemiology, and end results database. World J Oncol 2025; 16: 286.

31.

Bagaria

Heckman

Diehl

, et al. Delay to colectomy and survival for patients diagnosed with colon cancer. J Invest Surg 2019; 32: 350–357.

32.

Grass

Behm

Duchalais

, et al. Impact of delay to surgery on survival in stage i–iii colon cancer. Eur J Surg Oncol 2020; 46: 455–461.

33.

Kucejko

Holleran

Stein

, et al. How soon should patients with colon cancer undergo definitive resection?. Diseases Colon Rectum 2020; 63: 172–182.

34.

Stacy

. A generalization of the gamma distribution. Ann Math Stat 1962; 33: 1187–1192.

35.

Cain

Robins

Lanoy

, et al. When to start treatment? a systematic approach to the comparison of dynamic regimes using observational data. Int J Biostat 2010; 6: 18.

36.

Cole

Hernán

. Constructing inverse probability weights for marginal structural models. Am J Epidemiol 2008; 168: 656–664.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

2.16 MB

0.00 MB