Survival Analysis: Employee Attrition mit lifelines
Code
import numpy as npimport pandas as pdimport matplotlib.colors as mcolorsimport plotly.express as pximport plotly.graph_objects as goimport plotly.io as piopio.renderers.default ="notebook"from lifelines import KaplanMeierFitter, CoxPHFitterfrom lifelines.statistics import logrank_testpd.set_option("display.max_columns", 50)# Plot-Stil anpassen: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_templatepio.templates.default ="custom"DATA_PATH ="data/WA_Fn-UseC_-HR-Employee-Attrition.csv"
Dieses Notebook demonstriert Survival-Analyse-Methoden (Kaplan-Meier, Log-Rank-Test, Cox Proportional-Hazards) angewendet auf Mitarbeiterfluktuation (Attrition).Wir verwenden dazu einen öffentlichen Beispieldatensatz von Kaggle:IBM HR Analytics Employee Attrition & Performance
Der Datensatz
Die Kernidee dieser Analyse ist, dass wir eine Schätzung dazu haben wollen, wie lange Mitarbeitende im Unternehmen bleiben, abhängig von bestimmten Faktoren. Hauptproblem dabei sind rechtszensierte Daten: Für viele Individuen wissen wir garnicht, nach wie vielen Jahren sie kündigen werden, nur, dass sie mindestens von ihrer Einstellung bis heute dabei sind. Ein kleiner Auszug aus den Daten:
Bevor wir anfangen zu modellieren, schauen wir uns die Verteilung der duration (Betriebszugehörigkeit) und die Gesamt-Attritionsrate an. Wichtig: Personen mit event_observed = 0 sind rechtszensiert. Sie haben nicht gekündigt, sondern die Beobachtung endete, während sie noch angestellt waren. Sie sind mindestens bis zu dem Zeitpunkt geblieben, als wir gezählt haben, mehr können wir aber nicht über sie sagen. Sie einfach aus der Analyse auszuschließen oder als „werden nie kündigen“ zu behandeln, würde aber die Ergebnisse systematisch verzerren. Genau das vermeidet die Survival-Analyse.
Code
fig = ( px.histogram(df.duration) .update_layout( bargap=.1, xaxis_title="Jahre", yaxis_title="Anzahl Angestellte", title="<b>Verteilung der Verweildauer insgesamt</b><br>Wie viele Mitarbeitende sind oder waren X Jahre im Betrieb?", margin=dict(t=80, r=0, l=40), ) .update_traces( showlegend=False, hovertemplate="<b>Verweildauer</b><br>%{y} Mitarbeitende sind oder waren<br>%{x} Jahre bei uns.<extra></extra>", marker_color="#708ebb" ))fig.show()
Abbildung 1: Die Verteilung der Jahre-im-Betrieb aller Mitarbeitenden, ob ehemalige oder gegenwärtige.
Abbildung 2: Die gleiche Verteilung wie oben, aber nach Gruppen aufgeschlüsselt. Die beiden senkrechten Linien markieren jeweils den Mittelwert der Verweildauer über alle Mitglieder der jeweiligen Gruppe hinweg. Hier wird deutlich: Würden wir uns nur auf diejenigen konzentrieren, die wirklich gegangen sind, würden wir die durchschnittliche Verweildauer aller Mitarbeitenden unterschätzen.
Schätzung der mittleren Verweildauer
Wir sind zunächst daran interessiert, wie lange Mitarbeitende im Unternehmen bleiben. Genauer: Wie sich mit der Zeit die Wahrscheinlichkeit entwickelt, dass sie noch da sind. Die Zählungen oben legen eine Intuition nahe: Man könnte mit Blick auf die Verteilung der tatsächlich erfolgten Kündigungen (grün) über die Jahre versucht sein zu sagen: „Wenn wir wissen wollen, nach welcher Zeit eine 50%-Wahrscheinlichkeit besteht, dass jemand geht, dann suchen wir uns den Punkt auf der Zeitachse, der genau 50% der gekündigten Mitarbeitenden links und rechts von sich hat.“
Der Gedanke dahinter: Ohne zusätzliche Informationen gehen wir davon aus, dass die Zukunft nicht überraschend von der Vergangenheit abweicht und somit der Mittelwert der Verweildauer eine:r Mitarbeiter:in in der Vergangenheit unser best guess über die zu erwartende Verweildauer von Neuzugängen ist.
Diese Intuition hat ein Problem: Alljene Mitarbeitenden, die heute, zum Zeitpunkt der Zählung, aktiver Teil der Belegschaft sind. Erstens wissen wir prinzipiell nicht, wie lang ihr Verbleib sein wird. Und zweitens: Wie die beiden Linien für die Erwartunswerte der Verweildauer zeigen, würde uns eine Schätzung nur anhand der Ehemaligen in die Irre führen.
Kaplan-Meier-Schätzung
Der Kaplan-Meier-Schätzer ist ein statistisches Instrument, das den gesamten Datensatz inklusive der aktiven Mitarbeitenden berücksichtigt. Er gibt an, mit welcher Wahrscheinlichkeit ein beliebiges Individuum nach einer bestimmten Zeit t noch im Unternehmen ist, unter korrekter Berücksichtigung der Zensierung.
fig = go.Figure()add_km_ci(fig, km, "blue")fig.update_layout( hovermode="x", title="<b>Kaplan-Meier-Kurve:</b><br>Wahrscheinlichkeit, im Unternehmen zu bleiben", yaxis_title="Verbleibenswahrscheinlichkeit", xaxis_title="Jahre im Unternehmen", margin=dict(t=80, r=0, l=40), yaxis_range=[0, 1],)fig.show()
Abbildung 3: Die Kaplan-Meier-Kurve ist im Grunde eine alternative Darstellungsweise des Datensatzes mit zusätzlicher Bezifferung der Vorhersageunsicherheit. Das blaue Band zeigt ein 95%-iges Konfidenzintervall, das wiedergibt, wie groß die Unsicherheit über das Ergebnis ist. Nach hinten wird dieses Intervall breiter, weil hier immer weniger Daten vorliegen - nur wenige Mitarbeitende sind über 30 Jahre dabei.
Blick auf einzelne Gruppen
Mit der Analyse wie sie bis jetzt steht sehen wir den Schätzwert über die gesamte Belegschaft hinweg. Nun haben wir aber noch einige weitere Daten, die es uns erlauben, zu fragen, wie sich diese Kurve verhält, wenn wir die Belegschaft aufteilen in jene, die Überstunden auf dem Konto haben, gegenüber den anderen.
Code
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-Kurve</b><br>für die Gruppen ohne (orange) und mit (blau) Überstunden", yaxis_title="Verbleibenswahrscheinlichkeit", xaxis_title="Jahre im Unternehmen", yaxis_range=[0, 1], hovermode="x", margin=dict(t=80, r=0, l=40),)# Log-Rank-Test: Unterscheiden sich die beiden Kurven statistisch signifikant?group_yes = df[df["OverTime"] =="Yes"]group_no = df[df["OverTime"] =="No"]fig.show()
Abbildung 4: Zwei separate Kaplan-Meier-Kurven, nachdem wir die Belegschaft in gekündigt und noch aktiv aufgeteilt haben. Die Konfidenzintervalle sind breiter, weil pro Kurve weniger Daten verwendet werden. Die Bereiche, in denen sich die beiden Bänder nicht berühren zeigen, dass sich die Gruppen durchaus voneinander unterscheiden.
Das Ergebnis: Mitarbeitende mit einem aktiv gefüllten Überstundenkonto neigen eher dazu, früher aus dem Unternehmen zu gehen.
Code
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-Wert = {result.p_value:.4f}")if result.p_value <0.05:print("→ Statistisch signifikanter Unterschied zwischen den Gruppen (p < 0.05).")else:print("→ Kein statistisch signifikanter Unterschied zwischen den Gruppen (p >= 0.05).")
Log-Rank-Test (OverTime Yes vs. No): p-Wert = 0.0000
→ Statistisch signifikanter Unterschied zwischen den Gruppen (p < 0.05).
Und ein Log-Rank-Test drückt als Kennwert aus, was die Grafik bereits suggeriert.
Das Cox Proportional-Hazards-Modell
Die Kaplan-Meier-Kurven eignen sich gut, um den Effekt einer Variable auf einmal zu untersuchen. Wir teilen die Stichprobe in zwei Gruppen und schauen uns die jeweilige Kaplan-Meier-Kurve an. In der Realität wirken aber meist mehrere Faktoren gleichzeitig, und sie sind oft nicht unabhängig voneinander: Wer viel reist, verdient tendenziell auch mehr; wer unzufrieden ist, hat vielleicht schon öfter den Arbeitgeber gewechselt. Wenn wir immer nur einen einzigen Schnitt durch die Gesamtmenge machen können, können wir solche Effekte nicht auseinanderhalten.
Das Cox Proportional-Hazards-Modell schätzt dagegen den Effekt mehrerer Einflussfaktoren gleichzeitig. Dafür formuliert es die
Verbleibenswahrscheinlichkeit x Jahre nach Beginn
um in das
momentane Kündigungsrisiko zur Zeit t, nachdem jemand bis t geblieben ist.
Außerdem teilt es das Risiko auf in den Baseline-Hazard über die ganze Gruppe hinweg und jede multiplikative Verschiebung davon durch die Kovariaten (Überstunden, …). Wie groß die positive oder negative Verschiebung für eine bestimmte Kovariate ist, ermitteln wir durch das model fitting. Dabei wird das Verhältnis der Hazard-Raten zweier Individuen als über die Zeit konstant angenommen, daher der Name „proportional hazards“.
Auswahl und Aufbereitung der Kovariaten
Wir wählen eine Teilmenge an Merkmalen, die inhaltlich plausible Treiber von Fluktuation sind, und bereiten sie für das Modell auf:
MonthlyIncome ist stark rechtsschief verteilt, daher loggen wir es (logMonthlyIncome). So entspricht der geschätzte Koeffizient einer relativen statt einer absoluten Gehaltsänderung.
OverTime codieren wir als 0/1.
MaritalStatus und BusinessTravel sind kategorial ohne natürliche Ordnung und werden per One-Hot-Encoding in Dummy-Variablen überführt (jeweils eine Referenzkategorie wird weggelassen: Divorced bzw. Non-Travel).
Ordinalskalierte Zufriedenheits- und Balance-Werte (JobSatisfaction, EnvironmentSatisfaction, WorkLifeBalance, StockOptionLevel) gehen als numerische Stufen ein.
Wenn wir zwei zufällige Mitarbeitende auswählen und vorhersagen würden, wer zuerst geht, wäre dies ohne weite Daten ein Münzwurf mit 50 % Wahrscheinlichkeit. Der Concordance-Index von rund 0,86 sagt nun aus: Wenn wir das Modell dafür benutzen würden, anhand der Prädiktoren der beiden Individuen zu schätzen, dann lägen wir in 86 % der Fälle richtig. Das ist eine deutliche Verbesserung. Der Log-Likelihood-Ratio-Test gegen ein Nullmodell ohne Kovariaten ist mit p ≈ 2,5·10⁻⁸³ hochsignifikant: Die gewählten Merkmale erklären gemeinsam substanziell mehr als der Zufall.
Hazard Ratios im Überblick
Nun zu der großen Tabelle mit Informationen über jede einzelne Kovariate: Das sind eine Menge Zahlen, die sich besser visuell verstehen lassen und sich daher im Plot unten wiederfinden. Die Hazard Ratio auf der x-Achse lässt sich direkt interpretieren: Werte über 1 bedeuten ein erhöhtes, Werte unter 1 ein reduziertes Kündigungsrisiko pro Einheit der Kovariate, verglichen mit der Referenz. Der Forest Plot macht Effektstärke, Unsicherheit und Signifikanz auf einen Blick vergleichbar:
Abbildung 5: Die Hazard-Ratios der verschiedenen Kovariaten mit Effektstärke, Unsicherheit und farblich markierter Signifikanz, sortiert von stark positivem zu stark negativem Effekt auf den Hazard.
Die Ergebnisse decken sich mit der univariaten Analyse oben, quantifizieren sie aber und kontrollieren gleichzeitig für die jeweils anderen Merkmale:
Überstunden (OverTime): Hazard Ratio ≈ 3,22. Überstunden sind mit einem gut verdreifachten Kündigungsrisiko verbunden, der stärkste Einzeleffekt im Modell.
Einkommen (logMonthlyIncome): Eine Verdopplung des Monatsgehalts ist mit einem um rund 61 % reduzierten Risiko verbunden. Gehalt ist ein starker protektiver Faktor.
Reisetätigkeit: Gegenüber Non-Travel haben Vielreisende (Travel_Frequently) ein knapp verdreifachtes (HR ≈ 2,90), Gelegenheitsreisende (Travel_Rarely) ein rund 83 % höheres Risiko (HR ≈ 1,83, an der Signifikanzgrenze).
Familienstand: Singles kündigen deutlich häufiger als Geschiedene (Referenzgruppe, HR ≈ 2,12); Verheiratete unterscheiden sich statistisch nicht von der Referenz.
Zahl früherer Jobs (NumCompaniesWorked): +21 % Risiko pro vorherigem Arbeitgeber. Eine Historie häufiger Wechsel setzt sich fort.
Arbeitsweg (DistanceFromHome): leichter, aber signifikanter Effekt (+3 % pro Kilometer Arbeitsweg).
Alter (Age): -6 % Risiko pro Lebensjahr, ältere Mitarbeitende bleiben tendenziell länger.
Optionen der Firma (StockOptionLevel): Effekt in die erwartete (protektive) Richtung, aber nicht signifikant (p ≈ 0,22).
JobSatisfaction, EnvironmentSatisfaction, WorkLifeBalance: jede Stufe auf der 1–4-Skala reduziert das Risiko um jeweils rund 19–20 %, auch nachdem für Gehalt, Reisetätigkeit & Co. kontrolliert wurde.
Prüfung der Proportional-Hazards-Annahme
Das Modell setzt voraus, dass sich die Hazard-Raten zweier Ausprägungen einer Kovariate über die gesamte Zeitachse hinweg um einen konstanten Faktor unterscheiden. Ist das nicht der Fall, können Schätzungen verzerrt sein. lifelines bietet dafür einen statistischen Test (basierend auf skalierten Schoenfeld-Residuen):
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
Zwei Variablen (logMonthlyIncome und MaritalStatus_Married) fallen bei diesem Test knapp unter die 5-%-Schwelle, suggerieren also, dass ihr Einfluss auf die Verweildauer mit der Zeit linear variiert. Bei fast 1.500 Beobachtungen und einer zweistelligen Zahl an Kovariaten ist das nicht überraschend, schon aus reinem Zufall werden bei diesem Schwellenwert einzelne Merkmale auffällig, ohne dass die Modellschlüsse dadurch grundsätzlich infrage stehen. Für eine weiter gehende Produktivanalyse wäre der nächste Schritt, logMonthlyIncome zu binnen und darauf zu stratifizieren, oder einen Interaktionsterm mit der Zeit einzuführen. An dieser Stelle soll aber unsere Analyse enden.
Fazit
An den Kaplan-Meier-Kurven lässt sich ablesen, ob sich die Verbleibenswahrscheinlichkeit bestimmter Gruppen unterscheidet, etwa ob Überstunden einen Einfluss haben. Der Log-Rank-Test drückt zudem als Zahl aus, was der Plot der beiden Gruppen schon suggeriert: Hier gibt es tatsächlich einen Unterschied.
Das Cox-Modell bestätigt und schärft dieses Bild: Überstunden und geringe Zufriedenheit bleiben starke Treiber, auch wenn man gleichzeitig für Gehalt, Alter, Reisetätigkeit, Familienstand und Wechselhistorie kontrolliert. Gehalt erweist sich dabei als der stärkste protektive Hebel, dicht gefolgt von Zufriedenheits- und Work-Life-Balance-Werten, beides Stellschrauben, an denen ein Unternehmen aktiv drehen kann. Reisetätigkeit und ein langer Arbeitsweg sind eher strukturelle Risikofaktoren, die sich schwerer kurzfristig ändern lassen, aber bei der Gestaltung von Rollen und Standortentscheidungen mitgedacht werden sollten.
Aus dieser Analyse ließe sich ableiten, welche Stellhebel (z. B. Reduktion von Überstunden) das größte Potenzial haben, Fluktuation zu senken und für welche Mitarbeitendengruppen präventive Maßnahmen (Retention-Gespräche, Entwicklungsangebote) am wirkungsvollsten wären.
Damit ergänzen sich die beiden Methoden: Kaplan-Meier liefert die intuitiv verständliche, gruppenweise Visualisierung; Cox Proportional-Hazards liefern die multivariate, wechselseitig kontrollierte Quantifizierung. Beides zusammen ist die Grundlage für datengetriebene Retention-Strategien.