Piecewise Constant Hazards Calculations

Daniel Sabanés Bové


In this vignette we derive the explicit formula for the overall survival under piecewise constant hazards. This is implemented in the function PWCsurvOS().

So we are in the situation where the hazards for the transitions are piecewise constant, i.e. all three hazard functions \(\lambda_{01}(t)\) (stable to progression), \(\lambda_{02}(t)\) (stable to death) and \(\lambda_{12}(t)\) (progression to death) are step functions. Say the start points of the constant hazard pieces for \(\lambda_{01}(t)\) are \(0 \equiv t_{01}^{(1)} < \dotsb < t_{01}^{(k_{01})}\), \(k_{01} \geq 1\), with corresponding constant positive hazards \(h_{01}^{(1)}, \dotsc, h_{01}^{(k_{01})}\). Obviously we use here the smallest set of pieces, i.e. neighboring hazards are required to be different, \(h_{01}^{(j)} \neq h_{01}^{(j+1)}\). This holds analogously for the hazard functions of the other two state transitions.

We denote the cumulative hazards similarly as \(\Lambda_{01}(t)\), \(\Lambda_{02}(t)\) and \(\Lambda_{12}(t)\). Note that these are piecewise linear, with the slope changes occurring at the times of hazard changes.

Overall survival calculation

Now we want to calculate the overall survival (OS) survival function induced by the piecewise constant hazard model. We start from

\[ S_{\text{OS}}(t) = S_{\text{PFS}}(t) + \int_0^t S_{\text{PFS}}(u) \lambda_{01}(u)\exp(\Lambda_{12}(u) - \Lambda_{12}(t))\, du \] where \(S_{\text{OS}}(t)\) is the survival function for OS, and \(S_{\text{PFS}}(t)\) is the survival function for PFS with the closed form

\[ S_{\text{PFS}}(t) = \exp(- \Lambda_{01}(t) - \Lambda_{02}(t)). \] Hence we can rewrite the integral from above as

\[ \exp(- \Lambda_{12}(t)) \int_0^t \exp(\Lambda_{12}(u) - \Lambda_{01}(u) - \Lambda_{02}(u))\lambda_{01}(u)\, du \] So overall we now have

\[ S_{\text{OS}}(t) = S_{\text{PFS}}(t) + \exp(- \Lambda_{12}(t)) \cal{I}(t) \] and we can rewrite the integral \[ \cal{I}(t) := \int_0^t \exp(\Lambda_{12}(u) - \Lambda_{01}(u) - \Lambda_{02}(u))\lambda_{01}(u)\, du \] in terms of the unique starting time points \(0 \equiv t_{(1)} < \dotsb < t_{(k)}\), chosen such that the set \(\{t_{(1)}, \dotsc, t_{(k)}\}\) is the smallest super set of all state specific transition starting points \(\{t_{01}^{(1)}, \dotsc, t_{01}^{(k_{01})}\}\), \(\{t_{02}^{(1)}, \dotsc, t_{02}^{(k_{02})}\}\) and \(\{t_{12}^{(1)}, \dotsc, t_{12}^{(k_{12})}\}\), as

\[\begin{align} \cal{I}(t) = &\int_{t_{(1)}}^{t_{(2)}} \exp(a_{(1)} + b_{(1)}(u - t_{(1)}))h_{01(1)}\,du \\ &+ \dotsb + \\ &\int_{t_{(l)}}^{t} \exp(a_{(l)} + b_{(l)}(u - t_{(l)}))h_{01(l)}\,du, \end{align}\] where:

Note that this is essentially just because of \[ \Lambda(t) = \Lambda(s) + h (t-s) \] when there is a constant hazard \(\lambda(t) \equiv h\) and two time points \(s<t\).

We can then easily derive a closed form for each integral part, \(j = 1, \dotsc, l\), where \(t_{(l+1)}\equiv t\) is the end point of the last integral:

\[\begin{align} \cal{I}_{j} &= \int_{t_{(j)}}^{t_{(j+1)}} \exp(a_{(j)} + b_{(j)}(u - t_{(j)}))h_{01(j)}\,du\\ &= \exp(a_{(j)} - b_{(j)}t_{(j)})h_{01(j)} \int_{t_{(j)}}^{t_{(j+1)}} \exp(b_{(j)}u)\,du\\ &= \exp(a_{(j)} - b_{(j)}t_{(j)})h_{01(j)} \left[ b_{(j)}^{-1}\exp(b_{(j)}u) \right]_{t_{(j)}}^{t_{(j+1)}}\\ &= \frac{h_{01(j)}}{b_{(j)}} \exp(a_{(j)} - b_{(j)}t_{(j)}) (\exp(b_{(j)}t_{(j+1)}) - \exp(b_{(j)}t_{(j)}))\\ &= \frac{h_{01(j)}}{b_{(j)}} \exp(a_{(j)}) (\exp(b_{(j)}(t_{(j+1)} - t_{(j)})) - 1) \end{align}\]

Note that if it should happen that \(b_{(j)} = 0\), the integral simplifies further to \[ \cal{I}_{j} = \int_{t_{(j)}}^{t_{(j+1)}} \exp(a_{(j)})h_{01(j)}\,du = \exp(a_{(j)})h_{01(j)}(t_{(j+1)} - t_{(j)}). \]


The above formula is implemented in the function PWCsurvOS(). Note that there are a few modifications compared to above exposition:

  1. In order to be efficient, we look at a vector of time points \(t\) directly from the beginning. Then, find the unique time points and order them, and perform piecewise integration between the sorted unique time points.
  2. In order to avoid numerical overflow with large time points, we move the term \(\exp(- \Lambda_{12}(t))\) inside the integral, similarly as in the starting point of the derivation.