wrong result with CoxPHFitter

[takra]

Hello

I'm beginner in lifelines, so hope you will be indulgent. As a beginner I start to test Cox model in lifelines on known dataset of two groups of leukemia patients: one group of 21 persons has received a certain treatment; the other group of 21 persons has received a placebo. The data come from Freireich et al., Blood, 1963. Results of Cox model for this dataset is fully observed in David G. Kleinbaum, Mitchel Klein auth. Survival Analysis A Self-Learning Text.

In this book for single feature 'group' Cox coefficient is 1.509. But using lifelines I get coefficient 1.57 and I can't find reason for this. Besides, I solved the same problem in Matlab and in scikit-survival and got the same coefficient 1.509 as in book. Could you help me please to find my mistake in lifelines?

kleinbaum.zip

[CamDavidsonPilon]

Hi @takra, I downloaded your data and I got the same result in lifelines, 1.57. In R's survival, I also got 1.57. Can you post your Matlab and Python code for sci-kit survival? And also report what the final partial log-likelihood output is for those inferences?

For inference in such a simple dataset (small data, no covariates), things should be almost identical.

[takra]

I have corrected value from book and in matlab - 1.509 in first post. Besides I add matlab calculations and page from book. Look at model 1. For model 2 you will observe difference too.

kleinbaum.zip

[CamDavidsonPilon]

My guess is that the other methods are not converging "enough". For example, in lifelines and in R, a partial-log-lik or -85.01 is achieved. In the screenshot, a partial-log-lik of -86.38 is achieved. We want to maximize this value, so I suspect that the other methods are halting convergence too early.

I see the same situation happen in model 2: lifelines and R achieve a larger partial-log-lik, hence should be considered more accurate.

Can you change the convergence parameters in the Matlab code?

[takra]

Thank you for idea,
I have increased
coxphopt.MaxIter = 500;
coxphopt.TolX = 1e-9;
coxphopt.TolFun = 1e-9;
coxphopt.MaxFunEvals = 500 but loglikelihood doesn't change at all.

[CamDavidsonPilon]

hm hm hm. Those seem to be enough. It could be how ties are handled. Do you know if Matlab uses Breslow or Efron method for handling ties? I know that scikit survival doesn't use Breslow (which lifelines and R does) - that might be it!

EDIT: oops got that backwards: lifelines and R use Efron. scikit uses breslow

[CamDavidsonPilon]

yeaaa I just tested method="breslow"in R, and got the 1.5092 value.

So it's how ties are handled. Efron is considered the most accurate (vs performance), but take some time to read up and decide which tie-breaking method you want.

[takra]

By default matlab use breslow ties, after change to Efron get 1.57. Thank you very much!