Survival Analysis: Employee Attrition with lifelines

import numpy as np
import pandas as pd
import matplotlib.colors as mcolors
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = "notebook"

from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

pd.set_option("display.max_columns", 50)

# Adjust plot style:
custom_template = go.layout.Template(
    layout=go.Layout(
        paper_bgcolor="rgba(255, 255, 255, 0)",
        plot_bgcolor="#dddddd",
        title=dict(x=0.05, xanchor="auto"),
        yaxis_showline=False,
    )
)

pio.templates["custom"] = custom_template
pio.templates.default = "custom"

DATA_PATH = "data/WA_Fn-UseC_-HR-Employee-Attrition.csv"

This notebook demonstrates survival analysis methods (Kaplan-Meier, log-rank test, Cox proportional hazards) applied to employee attrition. We use a public example dataset from Kaggle for this: IBM HR Analytics Employee Attrition & Performance

The Dataset

The core idea of this analysis is that we want an estimate of how long employees stay with a company, depending on certain factors. The main challenge here is right-censored data: for many individuals we simply don’t know after how many years they will quit, only that they have stayed at least from their hire date until today. A small excerpt from the data:

# Derive target variables for survival analysis
df = pd.read_csv(DATA_PATH)
df["duration"] = df["YearsAtCompany"].clip(lower=0.01)  # duration > 0 required
df["event_observed"] = (df["Attrition"] == "Yes").astype(int)
df.head()

	Age	Attrition	BusinessTravel	DailyRate	Department	DistanceFromHome	Education	EducationField	EmployeeCount	EmployeeNumber	EnvironmentSatisfaction	Gender	HourlyRate	JobInvolvement	JobLevel	JobRole	JobSatisfaction	MaritalStatus	MonthlyIncome	MonthlyRate	NumCompaniesWorked	Over18	OverTime	PercentSalaryHike	PerformanceRating	RelationshipSatisfaction	StandardHours	StockOptionLevel	TotalWorkingYears	TrainingTimesLastYear	WorkLifeBalance	YearsAtCompany	YearsInCurrentRole	YearsSinceLastPromotion	YearsWithCurrManager	duration	event_observed
0	41	Yes	Travel_Rarely	1102	Sales	1	2	Life Sciences	1	1	2	Female	94	3	2	Sales Executive	4	Single	5993	19479	8	Y	Yes	11	3	1	80	0	8	0	1	6	4	0	5	6.00	1
1	49	No	Travel_Frequently	279	Research & Development	8	1	Life Sciences	1	2	3	Male	61	2	2	Research Scientist	2	Married	5130	24907	1	Y	No	23	4	4	80	1	10	3	3	10	7	1	7	10.00	0
2	37	Yes	Travel_Rarely	1373	Research & Development	2	2	Other	1	4	4	Male	92	2	1	Laboratory Technician	3	Single	2090	2396	6	Y	Yes	15	3	2	80	0	7	3	3	0	0	0	0	0.01	1
3	33	No	Travel_Frequently	1392	Research & Development	3	4	Life Sciences	1	5	4	Female	56	3	1	Research Scientist	3	Married	2909	23159	1	Y	Yes	11	3	3	80	0	8	3	3	8	7	3	0	8.00	0
4	27	No	Travel_Rarely	591	Research & Development	2	1	Medical	1	7	1	Male	40	3	1	Laboratory Technician	2	Married	3468	16632	9	Y	No	12	3	4	80	1	6	3	3	2	2	2	2	2.00	0

print(f"Number of observations: {len(df)}")
print(f"Observed events (resignations): {df['event_observed'].sum()} "
      f"({df['event_observed'].mean():.1%})")
print(f"Censored cases (still employed): {(1 - df['event_observed']).sum()} "
      f"({(1 - df['event_observed']).mean():.1%})")

Number of observations: 1470
Observed events (resignations): 237 (16.1%)
Censored cases (still employed): 1233 (83.9%)

Exploratory Analysis

Before we start modeling, let’s look at the distribution of duration (tenure) and the overall attrition rate. Important: individuals with event_observed = 0 are right-censored. They didn’t quit — the observation ended while they were still employed. They stayed at least until the point we counted, but we can’t say anything more about them. Simply excluding them from the analysis, or treating them as “will never quit,” would systematically bias the results. This is exactly what survival analysis avoids.

fig = (
    px.histogram(df.duration)
    .update_layout(
        bargap=.1,
        xaxis_title="Years",
        yaxis_title="Number of employees",
        title="<b>Overall distribution of tenure</b><br>How many employees have stayed or were staying for X years?",
        margin=dict(t=80, r=0, l=40),
    )
    .update_traces(
        showlegend=False,
        hovertemplate="<b>Tenure</b><br>%{y} employees have stayed or were staying<br>%{x} years with us.<extra></extra>",
        marker_color="#708ebb"
    )
)

fig.show()

Figure 1: The distribution of years-at-company across all employees, whether former or current.

df_grpdistro = df[["duration", "Attrition"]]

from plotly.subplots import make_subplots

fig = make_subplots(
    cols=1,
    rows=2,
    shared_yaxes=True,
    shared_xaxes=True,
    vertical_spacing=.01,
)

for atr, grp in df_grpdistro.groupby("Attrition"):
    for dur, dgrp in grp.groupby("duration"):

        fig.add_trace(
            go.Bar(
                x=dgrp.duration,
                y=[len(dgrp)],
                name=atr,
                marker_color={"Yes": "#759b75", "No": "#959595"}[atr],
                showlegend=False,
                marker_line_width=0,
                hovertemplate=f"<b>Tenure</b><br>%{{y}} employees have stayed or were staying<br>%{{x}} years with us.<extra>Event occurred: {atr}</extra>"
            ),
            col=1, row={"Yes": 1, "No": 2}[atr],
        )

    fig.add_trace(
        go.Scatter(
            x=[grp.duration.mean(), grp.duration.mean()],
            y=[0, 180],
            mode="lines",
            marker_color={"Yes": "#5d845d", "No": "#787878"}[atr],
            marker_line_dash="dot",
            showlegend=False,
            hoverinfo="skip",
        ),
        col=1, row={"Yes": 1, "No": 2}[atr],
    )

    fig.add_annotation(
        text="Former employees",
        x=.5, xanchor="center", xref="paper",
        y=170, yanchor="top", yref="y",
        showarrow=False,
        font_size=15,
    )

    fig.add_annotation(
        text="Currently employed (total tenure unknown)",
        x=.5, xanchor="center", xref="paper",
        y=170, yanchor="top", yref="y2",
        showarrow=False,
        font_size=15,
    )

fig.update_layout(
    bargap=.1,
    barmode="stack",
    title="<b>Distribution of tenure</b><br>split between former and currently employed.",
    yaxis_range=[0, 180],
    yaxis2_range=[0, 180],
    margin=dict(t=80, r=0, l=40),
    yaxis_title="Count", yaxis2_title="Count",
    xaxis2_title="Years"
)

fig.show()

Figure 2: The same distribution as above, but broken down by group. The two vertical lines each mark the mean tenure across all members of the respective group. This makes it clear: if we focused only on those who actually left, we would underestimate the average tenure of all employees.

Estimating Mean Tenure

We’re initially interested in how long employees stay with the company. More precisely: how the probability that they’re still there develops over time. The counts above suggest an intuition: looking at the distribution of resignations that actually occurred (green) over the years, one might be tempted to say: “If we want to know after how much time there’s a 50% probability that someone leaves, we just find the point on the time axis that has exactly 50% of the resigned employees to its left and right.”

The reasoning behind it: absent additional information, we assume the future won’t deviate surprisingly from the past, and so the mean tenure of employees in the past is our best guess for the expected tenure of new hires.

This intuition has a problem: all those employees who, at the time of counting, are an active part of the workforce. First, we fundamentally don’t know how long their tenure will be. And second: as the two lines for the expected-value tenures show, an estimate based only on former employees would mislead us.

Kaplan-Meier Estimation

The Kaplan-Meier estimator is a statistical tool that takes the entire dataset into account, including active employees. It states the probability that a given individual is still with the company after a certain time t, correctly accounting for censoring.

# Helper functions for repeated tasks

km = KaplanMeierFitter().fit(
    durations=df.duration,
    event_observed=df.event_observed,
    label="Total",
)

def to_rgba_str(color, alpha=0.3):
    r, g, b = mcolors.to_rgb(color)
    return f"rgba({int(r*255)}, {int(g*255)}, {int(b*255)}, {alpha})"

def add_km_ci(fig: go.Figure, km: KaplanMeierFitter, color: str, label: str = "Total") -> go.Figure:

    fillcolor = to_rgba_str(color)

    fig.add_trace(
        go.Scatter(
            x=km.timeline,
            y=km.confidence_interval_.iloc[:, 0],
            line_shape="hv",
            line_width=0,
            line_color=color,
            showlegend=False,
            legendgroup=label,
            name="Lower",
            hovertemplate="%{y:.2f}",
    ))

    fig.add_trace(
        go.Scatter(
            x=km.timeline,
            y=km.confidence_interval_.iloc[:, 1],
            line_shape="hv",
            line=dict(width=0, color=color),
            fill="tonexty",
            fillcolor=fillcolor,
            showlegend=False,
            legendgroup=label,
            name="Upper",
            hovertemplate="%{y:.2f}",
        ),
    )

    fig.add_trace(
        go.Scatter(
            x=km.timeline,
            y=km.survival_function_[label],
            mode="lines",
            line_shape="hv",
            line_color=color,
            showlegend=False,
            legendgroup=label,
            name=label,
            hovertemplate="%{y:.2f}",
        )
    )

    return(fig)

fig = go.Figure()
add_km_ci(fig, km, "blue")
fig.update_layout(
    hovermode="x",
    title="<b>Kaplan-Meier curve:</b><br>probability of staying with the company",
    yaxis_title="Probability of remaining",
    xaxis_title="Years with the company",
    margin=dict(t=80, r=0, l=40),
    yaxis_range=[0, 1],
)
fig.show()

Figure 3: The Kaplan-Meier curve is essentially an alternative way of presenting the dataset, with the added benefit of quantifying prediction uncertainty. The blue band shows a 95% confidence interval, reflecting how much uncertainty there is about the result. It widens toward the right because less and less data is available there — only few employees have been with the company for over 30 years.

A Look at Individual Groups

With the analysis as it stands so far, we see the estimate across the entire workforce. But we have some additional data that lets us ask how this curve behaves when we split the workforce into those who have logged overtime versus those who haven’t.

colors = {
    "Yes": "#228899",
    "No": "#ff8800",
}

fig = go.Figure()

for ovt, grp in df.groupby("OverTime"):
    km_group = KaplanMeierFitter()
    km_group.fit(
        grp["duration"],
        grp["event_observed"],
        label=f"OverTime = {ovt}"
    )
    add_km_ci(fig, km_group, colors[ovt], label=f"OverTime = {ovt}")

fig.update_layout(
    title="<b>Kaplan-Meier curve</b><br>for the groups without (orange) and with (blue) overtime",
    yaxis_title="Probability of remaining",
    xaxis_title="Years with the company",
    yaxis_range=[0, 1],
    hovermode="x",
    margin=dict(t=80, r=0, l=40),
)

# Log-rank test: do the two curves differ statistically significantly?
group_yes = df[df["OverTime"] == "Yes"]
group_no = df[df["OverTime"] == "No"]

fig.show()

Figure 4: Two separate Kaplan-Meier curves after splitting the workforce by overtime. The confidence intervals are wider because less data is used per curve. The regions where the two bands don’t touch show that the groups do in fact differ from one another.

The result: employees with actively accrued overtime tend to leave the company sooner.

result = logrank_test(
    group_yes["duration"], group_no["duration"],
    event_observed_A=group_yes["event_observed"], event_observed_B=group_no["event_observed"],
)
print(f"Log-rank test (OverTime Yes vs. No): p-value = {result.p_value:.4f}")
if result.p_value < 0.05:
    print("→ Statistically significant difference between the groups (p < 0.05).")
else:
    print("→ No statistically significant difference between the groups (p >= 0.05).")

Log-rank test (OverTime Yes vs. No): p-value = 0.0000
→ Statistically significant difference between the groups (p < 0.05).

And a log-rank test expresses as a single statistic what the chart already suggests.

The Cox Proportional-Hazards Model

Kaplan-Meier curves are well suited to examining the effect of one variable at a time. We split the sample into two groups and look at the respective Kaplan-Meier curve. In reality, however, several factors usually act at once, and they’re often not independent of one another: people who travel a lot also tend to earn more; people who are dissatisfied may have already changed employers more often. If we can only ever make a single cut through the whole dataset, we can’t disentangle such effects.

The Cox proportional-hazards model, by contrast, estimates the effect of several influencing factors simultaneously. To do so, it reformulates the

probability of remaining x years after start

into the

instantaneous risk of resignation at time t, given someone has stayed until t.

It also splits the risk into a baseline hazard shared across the whole group, and a multiplicative shift from that baseline for each covariate (overtime, …). How large the positive or negative shift is for a given covariate is determined through model fitting. The ratio of hazard rates between two individuals is assumed constant over time, hence the name “proportional hazards.”

Selecting and Preparing the Covariates

We select a subset of features that are plausible drivers of attrition on substantive grounds, and prepare them for the model:

• MonthlyIncome is strongly right-skewed, so we log it (logMonthlyIncome). This way the estimated coefficient corresponds to a relative rather than an absolute change in salary.
• OverTime is coded as 0/1.
• MaritalStatus and BusinessTravel are categorical with no natural order and are one-hot encoded into dummy variables (one reference category is dropped in each case: Divorced and Non-Travel respectively).
• Ordinal satisfaction and balance scores (JobSatisfaction, EnvironmentSatisfaction, WorkLifeBalance, StockOptionLevel) enter as numeric levels.

covariates = [
    "Age", "MonthlyIncome", "DistanceFromHome", "NumCompaniesWorked",
    "OverTime", "JobSatisfaction", "EnvironmentSatisfaction", "WorkLifeBalance",
    "StockOptionLevel", "MaritalStatus", "BusinessTravel",
]
df_cox = df[covariates + ["duration", "event_observed"]].copy()
df_cox["OverTime"] = (df_cox["OverTime"] == "Yes").astype(int)
df_cox["MonthlyIncome"] = np.log(df_cox["MonthlyIncome"])
df_cox = df_cox.rename(columns={"MonthlyIncome": "logMonthlyIncome"})
df_cox = pd.get_dummies(df_cox, columns=["MaritalStatus", "BusinessTravel"], drop_first=True)
bool_cols = df_cox.select_dtypes(bool).columns
df_cox[bool_cols] = df_cox[bool_cols].astype(int)

df_cox.head()

	Age	logMonthlyIncome	DistanceFromHome	NumCompaniesWorked	OverTime	JobSatisfaction	EnvironmentSatisfaction	WorkLifeBalance	StockOptionLevel	duration	event_observed	MaritalStatus_Married	MaritalStatus_Single	BusinessTravel_Travel_Frequently	BusinessTravel_Travel_Rarely
0	41	8.698347	1	8	1	4	2	1	0	6.00	1	0	1	0	1
1	49	8.542861	8	1	0	2	3	3	1	10.00	0	1	0	1	0
2	37	7.644919	2	6	1	3	4	3	0	0.01	1	0	1	0	1
3	33	7.975565	3	1	1	3	4	3	0	8.00	0	1	0	1	0
4	27	8.151333	2	9	0	2	1	3	1	2.00	0	1	0	0	1

Model Fit

First, some information from the model-fitting process:

cph = CoxPHFitter()
cph.fit(df_cox, duration_col="duration", event_col="event_observed")
cph.print_summary()

model	lifelines.CoxPHFitter
duration col	'duration'
event col	'event_observed'
baseline estimation	breslow
number of observations	1470
number of events observed	237
partial log-likelihood	-1360.05
time fit was run	2026-08-31 05:46:31 UTC

	coef	exp(coef)	se(coef)	coef lower 95%	coef upper 95%	exp(coef) lower 95%	exp(coef) upper 95%	cmp to	z	p	-log2(p)
Age	-0.07	0.94	0.01	-0.09	-0.04	0.92	0.96	0.00	-6.17	<0.005	30.43
logMonthlyIncome	-1.37	0.25	0.14	-1.65	-1.09	0.19	0.34	0.00	-9.59	<0.005	69.90
DistanceFromHome	0.03	1.03	0.01	0.02	0.05	1.02	1.05	0.00	4.20	<0.005	15.18
NumCompaniesWorked	0.19	1.21	0.03	0.14	0.24	1.15	1.27	0.00	7.42	<0.005	42.90
OverTime	1.17	3.22	0.13	0.91	1.43	2.48	4.18	0.00	8.82	<0.005	59.63
JobSatisfaction	-0.21	0.81	0.06	-0.33	-0.10	0.72	0.91	0.00	-3.61	<0.005	11.66
EnvironmentSatisfaction	-0.22	0.80	0.06	-0.33	-0.11	0.72	0.90	0.00	-3.79	<0.005	12.72
WorkLifeBalance	-0.22	0.80	0.09	-0.40	-0.04	0.67	0.96	0.00	-2.42	0.02	6.00
StockOptionLevel	-0.15	0.86	0.12	-0.39	0.09	0.68	1.10	0.00	-1.22	0.22	2.17
MaritalStatus_Married	0.23	1.26	0.22	-0.20	0.66	0.82	1.93	0.00	1.05	0.30	1.76
MaritalStatus_Single	0.75	2.12	0.27	0.23	1.28	1.26	3.58	0.00	2.82	<0.005	7.71
BusinessTravel_Travel_Frequently	1.06	2.90	0.32	0.44	1.69	1.56	5.40	0.00	3.35	<0.005	10.29
BusinessTravel_Travel_Rarely	0.61	1.83	0.30	0.01	1.20	1.01	3.32	0.00	2.00	0.05	4.47

Concordance	0.86
Partial AIC	2746.10
log-likelihood ratio test	428.15 on 13 df
-log2(p) of ll-ratio test	274.39

If we picked two random employees and tried to predict who would leave first, without any further data that would be a coin flip with a 50% probability. The concordance index of around 0.86 now tells us: if we used the model to estimate this from the two individuals’ predictors, we’d be right in 86% of cases. That’s a substantial improvement. The log-likelihood-ratio test against a null model without covariates is highly significant at p ≈ 2.5·10⁻⁸³: the chosen features jointly explain substantially more than chance.

Hazard Ratios at a Glance

Now for the large table of information about each individual covariate: that’s a lot of numbers, which are easier to grasp visually and are therefore shown in the plot below. The hazard ratio on the x-axis can be interpreted directly: values above 1 mean an increased, values below 1 a reduced risk of resignation per unit of the covariate, relative to the reference. The forest plot makes effect size, uncertainty, and significance comparable at a glance:

summary = cph.summary.sort_values("exp(coef)")

sig = summary["p"] < 0.05
colors = np.where(~sig, "#999999", np.where(summary["exp(coef)"] > 1, "#d7191c", "#1a9641"))

fig = go.Figure()

fig.add_vline(x=1, line_dash="dot", line_color="#666666")

fig.add_trace(go.Scatter(
    x=summary["exp(coef)"],
    y=summary.index,
    mode="markers",
    marker=dict(color=colors, size=10),
    error_x=dict(
        type="data",
        symmetric=False,
        array=summary["exp(coef) upper 95%"] - summary["exp(coef)"],
        arrayminus=summary["exp(coef)"] - summary["exp(coef) lower 95%"],
        color="#888888",
    ),
    hovertemplate="<b>%{y}</b><br>Hazard ratio: %{x:.2f}<extra></extra>",
))

fig.update_layout(
    title="<b>Cox model: hazard ratios</b><br>red = increased, green = reduced risk of resignation (grey = not significant, p ≥ 0.05)",
    xaxis_title="Hazard ratio (log scale)",
    xaxis_type="log",
    margin=dict(t=80, l=190, r=0),
)

fig.show()

Figure 5: The hazard ratios of the various covariates with effect size, uncertainty, and color-coded significance, sorted from strongly positive to strongly negative effect on the hazard.

The results align with the univariate analysis above, but quantify it while controlling for the other features at the same time:

• Overtime (OverTime): hazard ratio ≈ 3.22. Overtime is associated with a good tripling of resignation risk, the strongest single effect in the model.

• Income (logMonthlyIncome): doubling monthly salary is associated with a roughly 61% reduced risk. Salary is a strong protective factor.

• Business travel: compared to Non-Travel, frequent travelers (Travel_Frequently) have nearly triple the risk (HR ≈ 2.90), occasional travelers (Travel_Rarely) around 83% higher risk (HR ≈ 1.83, right at the significance threshold).

• Marital status: singles resign noticeably more often than divorced employees (reference group, HR ≈ 2.12); married employees don’t differ statistically from the reference.

• Number of previous jobs (NumCompaniesWorked): +21% risk per previous employer. A history of frequent changes continues.

• Commute (DistanceFromHome): a small but significant effect (+3% per kilometer of commute).

• Age (Age): -6% risk per year of age; older employees tend to stay longer.

• Company stock options (StockOptionLevel): effect in the expected (protective) direction, but not significant (p ≈ 0.22).

• JobSatisfaction, EnvironmentSatisfaction, WorkLifeBalance: each step on the 1–4 scale reduces risk by roughly 19–20%, even after controlling for salary, travel, and so on.

Checking the Proportional-Hazards Assumption

The model assumes that the hazard rates of two levels of a covariate differ by a constant factor across the entire time axis. If that’s not the case, estimates can be biased. lifelines provides a statistical test for this (based on scaled Schoenfeld residuals):

cph.check_assumptions(df_cox, p_value_threshold=0.05, show_plots=False);

The ``p_value_threshold`` is set at 0.05. Even under the null hypothesis of no violations, some
covariates will be below the threshold by chance. This is compounded when there are many covariates.
Similarly, when there are lots of observations, even minor deviances from the proportional hazard
assumption will be flagged.

With that in mind, it's best to use a combination of statistical tests and visual tests to determine
the most serious violations. Produce visual plots using ``check_assumptions(..., show_plots=True)``
and looking for non-constant lines. See link [A] below for a full example.

null_distribution	chi squared
degrees_of_freedom	1
model	<lifelines.CoxPHFitter: fitted with 1470 total...
test_name	proportional_hazard_test

		test_statistic	p	-log2(p)
Age	km	0.02	0.89	0.16
	rank	0.61	0.43	1.20
BusinessTravel_Travel_Frequently	km	0.22	0.64	0.65
	rank	0.20	0.65	0.62
BusinessTravel_Travel_Rarely	km	1.16	0.28	1.83
	rank	1.09	0.30	1.75
DistanceFromHome	km	0.31	0.58	0.80
	rank	0.05	0.82	0.28
EnvironmentSatisfaction	km	0.04	0.84	0.26
	rank	0.00	1.00	0.00
JobSatisfaction	km	1.53	0.22	2.21
	rank	0.66	0.42	1.26
MaritalStatus_Married	km	4.96	0.03	5.27
	rank	3.41	0.06	3.95
MaritalStatus_Single	km	1.09	0.30	1.75
	rank	0.72	0.40	1.33
NumCompaniesWorked	km	0.52	0.47	1.08
	rank	0.75	0.39	1.37
OverTime	km	0.31	0.58	0.79
	rank	0.00	0.98	0.03
StockOptionLevel	km	0.05	0.83	0.27
	rank	0.11	0.74	0.43
WorkLifeBalance	km	0.51	0.48	1.07
	rank	0.43	0.51	0.97
logMonthlyIncome	km	3.00	0.08	3.58
	rank	4.74	0.03	5.08

1. Variable 'logMonthlyIncome' failed the non-proportional test: p-value is 0.0295.

   Advice 1: the functional form of the variable 'logMonthlyIncome' might be incorrect. That is,
there may be non-linear terms missing. The proportional hazard test used is very sensitive to
incorrect functional forms. See documentation in link [D] below on how to specify a functional form.

   Advice 2: try binning the variable 'logMonthlyIncome' using pd.cut, and then specify it in
`strata=['logMonthlyIncome', ...]` in the call in `.fit`. See documentation in link [B] below.

   Advice 3: try adding an interaction term with your time variable. See documentation in link [C]
below.


2. Variable 'MaritalStatus_Married' failed the non-proportional test: p-value is 0.0259.

   Advice: with so few unique values (only 2), you can include `strata=['MaritalStatus_Married',
...]` in the call in `.fit`. See documentation in link [E] below.

---
[A]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html
[B]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Bin-variable-and-stratify-on-it
[C]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Introduce-time-varying-covariates
[D]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Modify-the-functional-form
[E]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Stratification

Two variables (logMonthlyIncome and MaritalStatus_Married) fall just below the 5% threshold in this test, suggesting their influence on tenure varies roughly linearly with time. With nearly 1,500 observations and a double-digit number of covariates, this isn’t surprising — by pure chance alone, individual features will stand out at this threshold without fundamentally calling the model’s conclusions into question. For a more thorough production analysis, the next step would be to bin logMonthlyIncome and stratify on it, or introduce an interaction term with time. At this point, however, our analysis will end here.

Conclusion

The Kaplan-Meier curves let us read off whether the probability of remaining differs between particular groups, for instance whether overtime has an effect. The log-rank test additionally expresses as a number what the plot of the two groups already suggests: there is in fact a difference here.

The Cox model confirms and sharpens this picture: overtime and low satisfaction remain strong drivers, even when controlling simultaneously for salary, age, travel, marital status, and job-change history. Salary turns out to be the strongest protective lever, closely followed by satisfaction and work-life balance scores — both levers a company can actively turn. Business travel and a long commute are more structural risk factors that are harder to change in the short term, but should be considered when designing roles and making location decisions.

From this analysis one could derive which levers (e.g., reducing overtime) have the greatest potential to lower attrition, and for which groups of employees preventive measures (retention conversations, development offers) would be most effective.

This is how the two methods complement each other: Kaplan-Meier provides the intuitively understandable, group-wise visualization; Cox proportional hazards provides the multivariate, mutually controlled quantification. Together, both form the basis for data-driven retention strategies.