Deriving survival from hazard function

Question

EDIT: It was a simple typo (+ instead of *), but the values still aren't quite right. Leaving the original in case others find it useful

I'm fairly new to working with survival functions so this might just simply be a misunderstanding of how the function works. I've been tasked with showing a manual calculation for Cox deriving the survival probability from the hazard function (so even though it exists, I can't just call the built in function). My understanding is:

The survival function should be derivable from the hazard function via

S(t)=exp{−∫h(x)dx} --from time = 0 to time = t

and the hazard function at any given time can be found using

h(t)=h0(t)*exp(b1x1+b2x2+...+bpxp) - where b is the coefficient for each covariate (first column)

I've used the lifelines package in Python to fit a cox-regression model to my data and gotten the covariates

    coef    exp(coef)   ...
DM  0.04    1.04    
HF  0.50    1.65    
PVD 0.29    1.34    
DEMENTIA    0.31    
AGE_AT_MOD_START    
MODALITY_NON-OPTIMAL START HD   -0.69   0.50    
MODALITY_OPTIMAL START HD   -1.36   0.26    
MODALITY_OPTIMAL START PD   -1.24   0.29    

And it has a baseline_hazard and baseline_cumulative_hazard function that can be called to get h0(t).

Lifelines has a built in function for calculating survival (predict_survival_function) which I've been using to check my work, but my answers have been quite off. I created an example patient and basically just plugged it into the h(t) equation above.

baseline hazard
0.0 0.000000
1.0 0.000000
2.0 0.000000
3.0 0.000000
4.0 0.000101
... ...
1822.0  0.077646
1823.0  0.077646
1824.0  0.077646
1825.0  0.000000
1826.0  0.000000

DM  HF  PVD DEMENTIA    AGE_AT_MOD_START    MODALITY_NON-OPTIMAL START HD   MODALITY_OPTIMAL START HD   MODALITY_OPTIMAL START PD
0   1   0   1   0       85                  1                              0                         0

h(4) = 0.000101+exp(0.04*0+0.50*1+0.29*0+0.31*1+0.04*85+....-1.24*0) = 20.90534 (uh oh)

I tried with cumulative baseline hazard as well, since baseline hazard is 0 for the days where no deaths occured, though my understanding is that the hazard formula just uses the baseline hazard, not the cumulative one. I'm positive the problem is inside the e^(covariates*variables) part- even just checking

h(0) = 20.90524

S(0) = 8.3362e-10

Is clearly wrong, since S(0) should be close to 1.0. The glaring potential offender is that I just plugged in age (85) directly since it was already numerical, but should all values be normalized or something? I've looked all over but I haven't seen an example of anyone doing the calculations, it's all just 'call a function' nowadays (which, fair, I'd do the same if I was allowed). Any insight would be appreciated.

For comparison, the build in survival function results

4.0 0.999724
6.0 0.999448
8.0 0.998897
9.0 0.998344
10.0    0.997792
... ...
1820.0  0.155034
1821.0  0.155034
1822.0  0.125425
1825.0  0.125425
1826.0  0.125425

Answer 1

I'm the author of lifelines. For prediction, we normalize the inputted values against the training dataset's mean.

So

h(4) = 0.000101 * exp(0.04*(0 - mean(DM))+0.50*(1 - mean(HF))+ ... ) 

Let me know if that clears things up.

---

Edit: This might be confusing for future readers, so let me make something more clear. We normalize by the mean interally, and this is abstracted from the end user. We do it for numerical stability reasons, but we "add it back" later. This doesn't not change the final prediction results.

Answer 2

The math on deriving the survival function from the hazard function is not that complex and I would recommend an exploration in a spreadsheet setting.

First, per this source, for example, confirming your equation, quoting:

${S_Y(y) = exp(−H_Y(y))}$

where:

${H_Y(y) = \int_0^y h_Y(t)dt}$

The function ${H_Y(y)}$ is called the cumulative hazard function or the integrated hazard function. Like the hazard function, the cumulative hazard function is not a probability. However, it is also a measure of risk: the greater the value of ${H_Y(y)}$, the greater the risk of failure by time ${y}$.

Now, we can approximate the cited integration by the use of numerical methods to enable the estimation of the survival function. Per this source, for example, to quote:

Trapezoidal rule of integration

The trapezoidal rule of integration approximates the function {f(x)} to be integrated by a first order polynomial {f_1(x) = a_0 + a_1 x}. If we chose {(a, f(a))} and {(b, f(b))} as the two points to approximate {f(x)} by {f_1(x)} we can solve for {a_0} and {a_1}, leading to

\begin{array}{rcl} I & \approx & \int_{a}^{b} (a_0 + a_1x)dx\\ & = & a_0(b-a) + a_1 \bigg(\frac{b^2 - a^2}{2}\bigg)\\ & = & \frac{(b-a)}{2}(f(a) + f(b)) \end{array}

To be more accurate, we usually divide the interval {[a, b]} into {n} sub-intervals, such that the width of each sub-interval is {h = (b - a)/n}. In this case, the trapezoidal rule of integration is given by

${\text {I } \approx \frac{(b-a)}{2n}[f(a) + 2(\sum_{i=1}^{n-1} f(a+ih)) + f(b)]}$

In the current case, you have empirical measures or estimates of the function, ${f(a+ih)}$ (or, per above ${h_Y(t)}$ ).