# Piecewise Constant Hazards Calculations

#### Daniel Sabanés Bové

Source:`vignettes/pwc_survival.Rmd`

`pwc_survival.Rmd`

## Introduction

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:

- \(t_{(l)}\) is the start of the last integral part and defined as the maximum starting point that is smaller than \(t\)
- \(h_{01(j)} = \lambda_{01}(t_{(j)})\), \(j=1, \dotsc, k\) are the hazard values for the stable to progression transition within the unique time intervals
- \(a_{(j)} = \Lambda_{12}(t_{(j)}) - \Lambda_{01}(t_{(j)}) - \Lambda_{02}(t_{(j)})\), \(j=1, \dotsc, k\) are the intercepts
- \(b_{(j)} = h_{12(j)} - h_{01(j)} - h_{02(j)}\), \(j=1, \dotsc, k\) are the slopes, again using the hazard values within the unique time intervals for the specific transitions

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)}). \]

## Implementation

The above formula is implemented in the function
`PWCsurvOS()`

. Note that there are a few modifications
compared to above exposition:

- 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.
- 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.