Add tutorial on how to estimate a Cox Proportional Hazards Model in glum

[MatthiasSchmidtblaicherQC]

And: Add link to Matt Mills' tutorial on fitting penalized splines in glum.

Checklist

• Added a CHANGELOG.rst entry

[stanmart]

Suggested change

"In the Cox model, the rate of event occurrence, $\\lambda(t,x_i)$, factorizes nicely into a linear predictor $\\eta_i=\\sum_k \\beta_k x_{ik}$ that depends on individual $i$'s characteristics but not on time $t$ and a baseline hazard $\\lambda_0$ that depends only on time, $\\lambda(t,x_i)=\\lambda_0(t)\\exp(\\eta_i)$ (the proportional hazards assumption). The partial log-likelihood of $\\eta_i$ is\n",

"In the Cox model, the rate of event occurrence, $\\lambda(t,x_i)$, factorizes nicely into a linear predictor $\\eta_i=\\sum_k \\beta_k x_{ik}$ that depends on individual $i$'s characteristics but not on time $t$, and a baseline hazard $\\lambda_0$ that depends only on time: $\\lambda(t,x_i)=\\lambda_0(t)\\exp(\\eta_i)$. This is known as the proportional hazards assumption). The partial log-likelihood of $\\eta_i$ is\n",

[stanmart]

Suggested change

"We now show that a Poisson approach in `glum` yields the same parameter estimates as a Cox model. For the latter, we use the [lifelines](https://github.com/CamDavidsonPilon/lifelines) library. We also take the dataset from lifelines, which is from an RCT on recidivism for 432 convicts released from Maryland state prisons with first arrest after release as event. We first load imports and the dataset. The dataset has one row per convict, with two outcome columns, the `week` until which the observation lasts and `arrest`, which indicates whether an arrest event happened or not (censoring)."

"We now show that a Poisson approach in `glum` yields the same parameter estimates as a Cox model. For the latter, we use the [lifelines](https://github.com/CamDavidsonPilon/lifelines) library. We also take the dataset from lifelines, which is from an RCT on recidivism for 432 convicts released from Maryland state prisons with first arrest after release as event. We first load imports and the dataset. The dataset has one row per convict, with two outcome columns: the `week` until which the observation lasts and `arrest`, which indicates whether an arrest event happened or not (censoring)."

[stanmart]

Suggested change

	"One might, therefore, wonder if the Poisson approach is competitive in terms of estimation speed. For the dataset here, the Poisson approach, including the data transformation by `survival_split`, turns out to be faster than the Cox model. This speedup is likely aided by tabmat's optimizations for the high-dimensional `week` categorical."
	"One might, therefore, wonder if the Poisson approach is competitive in terms of estimation speed. For the dataset here, the Poisson approach, including the data transformation by `survival_split`, turns out to be faster than the Cox model. This speedup is aided by tabmat's optimizations for the high-dimensional `week` categorical."

I'm confident this is the case :)

[stanmart]

Maybe mention here that there are various ways to deal with this? (e.g., exact - and expensive - calculations, efron approximation (reasonably accurate), breslow approximation (fastest). Also, in the case of many ties, an inherently discrete survival model might be the best option.

[stanmart]

They do use the same dataset-transformation (which they call stacking) idea, but they rely on a binomial model instead of Poisson and say that it is a good enough approximation. So relevant, but I would not say it's the same thing.

[stanmart]

Thank you Matthias, I like this tutorial a lot!

Maybe we could also mention that the model we estimate is conceptually a piecewise constant hazard model, but because we have a different constant for each event time, it becomes essentially non-parametric. As a corollary, if one wants to reduce the number of parameters, they can opt for merging some categories and estimating fewer "constants".