Survival Analysis and the Cox Proportional Hazards Model: A Beginner-Friendly Guide

ordinary regression can’t answer

Here is a dataset. 432 people were released from prison and were followed for one year. They were being watched for a single event: re-arrest. Some were re-arrested in week 8, some in week 30. Most (318 of them) reached the end of the study having never been re-arrested at all.

The simple question that was being asked in this study was how long until someone re-offends?

We cannot just average the arrest times. Three-quarters of the people never got arrested, so they have no arrest time to average. We cannot throw those people away either. Because if we only analyse the ones who failed, we will end up concluding that everyone eventually re-offends, which is both false and grim. And we also cannot code the survivors as arrested after week 52, because they weren’t. They were either arrested at some week after the 52nd week, that we never got to see, or maybe never.

That last situation, which is we know someone lasted at least this long, but we do not know how long in total, is called censoring. It is the entire reason survival analysis exists as its own field. Ordinary linear regression has no way to represent “the answer is at least 52.” It wants a number. Survival analysis is the set of tools built for the case where, for a big chunk of our data, the clock was still running when we stopped watching.

This post builds on that single problem. We will start with three core ideas we need to understand Survival Analysis. Then we will move on to the following three activities.

1. First, we will estimate a survival curve straight from data with Kaplan-Meier.

2. Then spend most of our time on the Cox proportional hazards regression.

3. And finally, we will fit one model in Python, reading its output as hazard ratios.

Three Core Concepts in Survival Analysis

Almost everything in survival analysis is built from three objects as detailed below.

1. The event and the time: Let’s pick one well-defined event. It can be death, re-arrest, machine failure, subscription cancellation or a loan default. For each subject, we record two things: how long they were observed, i.e., the duration, and whether the observation ended because the event happened or because we stopped watching, i.e., the event indicator, which is binary coded as 1 or 0. These two columns are the key information of every model below.

2. The survival function, S(t): This is the probability that a subject makes it past time t without the event. It starts at 1 (at time zero) and decays toward 0 over time. For example,

Survival at 12 months is 0.7″ means 70% are expected to still be event-free at a year.

3. The hazard function, h(t): The hazard is the instantaneous rate of the event at time t, given that we have survived up to t. Informally, we can say it like this:

Of the people who’ve made it this far, what fraction fail right now?

The distinction between S(t) and h(t) matters. Survival is cumulative; hazard is momentary. Our overall probability of still being alive at 80 is low (survival is small), but our hazard at the exact instant we turn 80 is a different quantity entirely. Survival is how much water is left in the tank. Hazard, on the other hand, is how fast it is draining right now. They are linked by a bit of calculus. Hazard is failure density divided by survival. So, integrating the hazard gives us back survival.

Now, you may think, why do we need to obsess over the hazard rather than just modelling survival directly? Because the hazard is the natural place to hang covariates. It’s easy to say “financial aid multiplies our rate of re-arrest by 0.68 at every moment.” That sentence is a statement about the hazard, and it’s exactly what the Cox model will give us.

Kaplan-Meier: a survival curve with no assumptions

Before we model anything, we can estimate S(t) directly from data. The Kaplan-Meier estimator (Kaplan & Meier, 1958, one of the most-cited papers in Statistics) does this without assuming any particular shape of the curve.

The idea is very simple. Let’s walk forward in time. At each moment where an event actually happens, we look at how many people were still at risk just before and how many failed. Multiply together those “fraction who survived this instant” numbers as we go. Censored people contribute to the at-risk count right up until they leave, and then quietly drop out without ever causing a downward step.

Now let’s apply this to our recidivism dataset. The specific dataset used in the article is the Rossi recidivism dataset, which is included with the lifelines Python library. It comes from a 1980 randomized experiment by Rossi, Berk, and Lenihan. The study followed 432 released prisoners for one year and recorded whether they were re-arrested, plus covariates like race, education etc. We will split the people into two groups: those who received financial aid after release, and those who did not. In the original study, this aid was assigned randomly, so the comparison is fair.

from lifelines import KaplanMeierFitter
from lifelines.datasets import load_rossi
import matplotlib.pyplot as plt

df = load_rossi()          # 432 rows: 'week', 'arrest' (1=event), + covariates
kmf = KaplanMeierFitter()

for value, label in [(0, "No financial aid"), (1, "Financial aid")]:
    g = df[df.fin == value]
    kmf.fit(g["week"], g["arrest"], label=label)
    kmf.plot_survival_function()
plt.ylabel("S(t): probability of remaining un-arrested")
plt.show()

Figure 1. Kaplan-Meier curves by financial-aid group. Every downward step is a re-arrest; the shaded bands are 95% confidence intervals. The aid group (blue) stays higher throughout. By week 52, about 22% of the aid group had been re-arrested versus 31% of the no-aid group.

To test whether that gap is more than noise, the standard tool is the log-rank test. It compares the observed number of events in each group against what we would expect if the two curves were really the same.

from lifelines.statistics import logrank_test

results = logrank_test(
    df[df["fin"] == 1]["week"],
    df[df["fin"] == 0]["week"],
    event_observed_A=df[df["fin"] == 1]["arrest"],
    event_observed_B=df[df["fin"] == 0]["arrest"],
)
results.p_value

On this data, it returns p ≈ 0.05, which is right on the border. This is a good reminder that Kaplan-Meier plus log-rank is a description of one variable at a time. It can not adjust for age, prior record, or anything else. The moment you want to control for covariates, we need a Cox model.

The Cox model: Regression on the hazard

When we want something like regression, e.g., plug in covariates, get out their effects, but for the hazard, the naive move is to write down a full formula for h(t) and estimate everything. But h(t) has a shape over time, and we usually have no idea what that shape is and no real desire to commit to one. Cox’s insight was that we can estimate the covariate effects without ever specifying the baseline shape. The model is specified as shown below:

Read it as two pieces multiplied together:

• h₀(t), the baseline hazard: the risk over time for a hypothetical subject with all covariates at zero. The shape is left completely unspecified. It can be any shape. This is the “nonparametric” part.
• exp(βᵀx), the covariate effect: a single number that scales the entire baseline up or down depending on personal characteristics. This is the “parametric” part.

Because the model is part-nonparametric and part-parametric, it is called a semiparametric model.

Here is the part that makes Cox model very useful – We can take two subjects and form the ratio of their hazards. The baseline h₀(t) is identical for both, so when we form the ratio, that baseline appears on the top and bottom and cancels out:

The right-hand side has no t in it. The unknown, time-varying baseline is gone. So, now the hazard ratio is the same at week 1, week 20, and week 52. We never had to specify the baseline’s shape. That cancellation is the entire magic trick.

Cox then turned this into a method called the partial likelihood. At each event time, the model asks a question:

Among everyone still at risk right now, how much more likely was it that this particular person failed, rather than one of the others?

Censored people fit naturally into this setup. They stay in the at-risk set until they leave the study, then quietly drop out. They never need an event. They still contribute information by telling us they were event-free up to that point.

Two things worth naming before we fit it:

• Ties : The partial likelihood assumes we can order the events. When two events land in the same week, software applies a correction (usually Efron’s method, which is the modern default and more accurate than Breslow’s).
• exp(β) is the hazard ratio. A coefficient of β = 0 means exp(β) = 1, i.e. no effect. β < 0 means exp(β) < 1, a protective factor that lowers the hazard. β > 0 raises it. We almost always report the exponentiated version.

Fitting the Cox model in Python and interpreting hazard ratios

The lifelines library makes fitting a Cox model straightforward.

from lifelines import CoxPHFitter

cph = CoxPHFitter()
cph.fit(df, duration_col="week", event_col="arrest")
cph.print_summary()

The estimated hazard ratios are presented below:

Covariate	Hazard ratio exp(β)	95% CI	p-value
Financial aid	0.68	0.47 – 1.00	0.047
Age (per year)	0.94	0.90 – 0.99	0.009
Prior convictions (each)	1.10	1.04 – 1.16	0.001
Race	1.37	0.75 – 2.50	0.31
Work experience	0.86	0.57 – 1.31	0.48
Married	0.65	0.31 – 1.37	0.26
On parole	0.92	0.63 – 1.35	0.67

To summarize:

• Financial aid, HR 0.68. Receiving aid is associated with a 32% lower hazard of re-arrest at any given moment. This is the treatment effect the original experiment was built to find.
• Age, HR 0.94. Each additional year of age lowers the hazard by about 6%. Older releasees re-offend more slowly.
• Prior convictions, HR 1.10. Each prior raises the hazard about 10%, and it multiplies: five priors is roughly 1.10⁵ ≈ 1.6× the hazard.

Hazard ratio is not a difference in survival probability. It is a multiplier on the momentary rate, assumed constant across the whole follow-up. This raises our next question: is it actually constant? The rest of the post explains this.

The assumption in the name: Proportional hazards

The model is called proportional hazards for a reason. Let’s look again at the ratio that made everything work:

The right-hand side has no ‘t’ in it. That is the proportional hazards (PH) assumption stated precisely: the hazard ratio between any two subjects is constant over time. Financial aid lowers the hazard by 32% in week 2 and by exactly 32% in week 50. The two groups’ hazards move up and down together. One is always a fixed multiple of the other. Their survival curves can never cross.

Sometimes that is true. Often it is not. A treatment might help enormously at first and wear off. A risk factor might only bite in the long run. When the real hazard ratio drifts with time, a plain Cox model averages it into a single number that’s wrong at both ends. So, we have to check.

The check most people should use is based on Schoenfeld residuals (Schoenfeld, 1982; and the scaled version with its formal test from Grambsch & Therneau, 1994). At each event time, the Schoenfeld residual measures the gap between the covariate value of the person who actually failed and the average covariate value among everyone at risk. If the PH assumption holds, those residuals should show no trend against time. If we can see a slope, the effect is changing over time, and proportionality is broken.

cph.check_assumptions(df, p_value_threshold=0.05, show_plots=True)

On our data the test flags two clear violations. Variable ‘age’ and ‘wexp’ failed the non-proportional test with p-values 0.0007 and 0.0063 respectively.

This example shows the real result on a real, famous dataset. It is exactly the kind of thing that stays invisible if we fit the model, read the p-values, and walk away. The coefficient for age (HR 0.94) isn’t wrong but incomplete. It is basically collapsing a genuinely time-varying effect into one number.

What to do when proportionality breaks

Finding a violation isn’t a dead end. It is usually the most interesting finding in the analysis, and there are three standard remediation strategies.

1. Stratify: If a variable violates PH but we only need to adjust for it (not estimate its effect), then we can put it in a strata. Stratification fits a separate baseline hazard for each level of that variable and never forces its effect to be proportional.

cph_strat = CoxPHFitter()
cph_strat.fit(
    df,
    duration_col="week",
    event_col="arrest",
    strata=["wexp"]
)
cph_strat.print_summary()

After stratifying on work experience to fix the proportional hazards violation, the key results hold i.e., financial aid (HR 0.68), older age (HR 0.94 per year), and prior convictions (HR 1.09 each) remain statistically significant, while race, marital status, and parole status are not. The model’s concordance of 0.61 indicates modest ability to rank who gets re-arrested sooner.

2. Let the effect vary with time. If we actually care how the effect changes, we need to add an interaction between the covariate and a function of time. This lets its hazard ratio rise or fall over follow-up. In lifelines this is CoxTimeVaryingFitter with the data split at event times. Basically, we are fitting β(t) instead of a single β.

3. Check the functional form first. The Schoenfeld test is sensitive to a misspecified covariate. If age’s true effect is nonlinear (risk falling fast then leveling off), entering it as a straight line can trip the PH test even when hazards are genuinely proportional. Before reaching for time-varying models, try a squared term on the offending variable.

Key takeaways

Censoring is a feature, not missing data. The whole point of these methods is to use the partial information in censored subjects correctly. If you ever find yourself dropping the people who didn’t have the event, you are basically introducing exactly the bias the field was invented to avoid.

Kaplan-Meier to see, Cox to adjust. Plot the KM curves first. They are assumption-light and they build intuition for the shape of survival. Move to Cox when you need to control for covariates. Report hazard ratios with confidence intervals, not just p-values.

A hazard ratio multiplies your risk, it doesn’t subtract from it. An HR of 0.68 means aid cuts the rate of re-arrest to about two-thirds at every moment. It doesn’t tell you how many people avoid arrest, or how much longer they stay free. Also, it quietly assumes that two-thirds holds the whole time, which is the thing you have to check.

Always test proportional hazards. It is two lines of code and it is the single thing that most separates a defensible survival analysis from a fragile one. When it fails, you can stratify, model the time variation, or fix the functional form.

Concordance, not R², for fit. The concordance index measures how well the model ranks who fails sooner.

Survival analysis is really just regression that’s honest about what it doesn’t know yet. Master the hazard ratio, respect the assumption in its name, and you’ll get answers that hold up long after the study ends.

References

Foundational

• Cox, D. R. (1972). Regression Models and Life-Tables. Journal of the Royal Statistical Society, Series B, 34(2), 187–220.
• Kaplan, E. L., & Meier, P. (1958). Nonparametric Estimation from Incomplete Observations. Journal of the American Statistical Association, 53(282), 457–481.

Diagnostics and the Cox model in depth

• Therneau, T. M., & Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer.
• Grambsch, P. M., & Therneau, T. M. (1994). Proportional Hazards Tests and Diagnostics Based on Weighted Residuals. Biometrika, 81(3), 515–526.
• Schoenfeld, D. (1982). Partial Residuals for the Proportional Hazards Regression Model. Biometrika, 69(1), 239–241.