0
$\begingroup$

I'm using the Cox regression model in lifelines (python) to try and predict what the probability of a patient surviving X days is given several variables.

Do to some very silly restrictions on the program a team is using, they need all the information in a flat file, like a .csv, so I can't use the built in functions of the package. I have to do the math manually. My initial thought was, if I just save all the coefficients (including baseline hazard for the timeframe they're interested in) I should be able to just plug it into the formula and get an output.

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

After fitting the model, I got the coefficients for each of my variables

coef    exp(coef)   se(coef)    coef lower 95%  coef upper 95%  exp(coef) lower 95% exp(coef) upper 95% z   p   -log2(p)
AGE_AT_MOD_START    0.04    1.04    0.00    0.03    0.04    1.03    1.04    17.68   <0.005  229.94
MODALITY_OPTIMAL START HD   -0.66   0.51    0.06    -0.77   -0.55   0.46    0.58    -11.75  <0.005  103.46
MODALITY_OPTIMAL START PD   -0.54   0.58    0.08    -0.69   -0.40   0.50    0.67    -7.19   <0.005  40.48
MODALITY_PLANNED NO DIALYSIS    0.69    1.99    0.08    0.54    0.84    1.72    2.31    9.10    <0.005  63.24
DM_Y    0.03    1.03    0.06    -0.09   0.16    0.91    1.17    0.50    0.62    0.70
HF_Y    0.50    1.65    0.05    0.40    0.60    1.49    1.82    9.92    <0.005  74.61
PVD_Y   0.30    1.35    0.05    0.19    0.40    1.21    1.49    5.67    <0.005  26.09
DEMENTIA_Y  0.30    1.35    0.09    0.12    0.47    1.13    1.60    3.37    <0.005  10.35
Concordance 0.75
Log-likelihood ratio test   1403.50 on 8 df
-log2(p) of ll-ratio test   986.63

But it fell apart trying to do the math manually. I wanted to give it a try to see what the probability of a patient dying on day 365 (I know to get survival I'll need to integrate from 0 to 365, but I wanted to start with the simpler stuff)

h(365) = h0(365) + exp(...coeff*covariates...)

But when I try to get the base hazard at that time

cph.baseline_hazard_.loc[365, :-1]

I get h0(365) = 0.0. Looking at the values before and after it

cph.baseline_hazard_.loc[360:370, :]

baseline hazard
360.0   0.000300
361.0   0.000150
362.0   0.000300
363.0   0.000301
364.0   0.000301
365.0   0.000000
366.0   0.000452
367.0   0.000454
368.0   0.000607
369.0   0.000152
370.0   0.000304

It's clear that there is a base hazard. If it's zero, should I just use the last nonzero value before it in my calculations? What the zero itself represents confuses me a little too, since it's obvious that it shouldn't mean that a patient's chance of dying is 0% on day 365.

If I use lifelines given survival function prediction and give it some patient's data and plot it, I can see that at day 365 there is definitely a number

tr_rows = cox_dummy.drop(['CVD_DAYS','DEADORALIVE'], axis=1)
tr_rows
cph.predict_survival_function(tr_rows.head()).plot() 

Cox survival plot

This is my first time working with any sort of survival analysis, so apologies if I'm missing something obvious.

$\endgroup$

2 Answers 2

2
$\begingroup$

EdM brings up some good points, please read them carefully. The solution to your problem is to not use the baseline_hazard_, but the cumulative_baseline_hazard_

cph.cumulative_baseline_hazard_.loc[365, :-1] * np.exp(np.dot(X, cph.params))

in some R survival programs it represents a set of some "average" values of the predictors. So read the manual carefully.

This is what lifelines does, so you need to "demean" the matrix X first.

$\endgroup$
1
$\begingroup$

The baseline hazard function in a Cox semi-parametric survival regression is empirical, with 0 hazard between observed event times. If no one died on day 365 in this cohort, then empirically there was 0 probability of death on that date.

What's shown in your plot is the the survival curve representing the integrated effect of the hazard over time. Can't see much detail right around day 365, and it's possible that cph.predict_survival_function smooths its output, too (don't use Python for this myself). The empirical survival curve, at least, would be flat between days 364 and 366.

So I fear that your attempt to make things "simpler" and your specific choice of day 365 is leading to some confusion instead in this case.

Other warnings:

Be careful about just how the baseline hazard is defined in your Python program. It might not be the hazard at reference values of all predictors; in some R survival programs it represents a set of some "average" values of the predictors. So read the manual carefully.

With the data that you show you will not get reliable confidence intervals for your estimates. As with any multiple regression model, you need the full covariance matrix of coefficient estimates, not just the standard errors of the individual coefficients, to get confidence intervals. See for example Wikipedia on the variance of a sum of correlated variables, as the regression coefficients estimates $\beta_i$ involved in the sum to get the linear predictor ($\sum \beta_i x_i$) are almost certainly not uncorrelated.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.