Review-Tutorial Survival Analysis

Lernziele

Nach dem Erarbeiten des Tutorials sind Sie in der Lage …

die grundlegenden Konzepte der Survival Analysis (Überlebenszeit, Überlebensfunktion, Hazard-Rate und Zensierung) zu erläutern.
Überlebensdaten in R aufzubereiten und mit einer Zeitvariable sowie einem Ereignisindikator korrekt abzubilden.
den Kaplan-Meier-Schätzer aus bedingten Überlebenswahrscheinlichkeiten zu erklären, in R anzuwenden und die geschätzten Überlebenskurven zu interpretieren.
die Idee des Log-Rank-Tests aus dem Vergleich beobachteter und unter der Nullhypothese erwarteter Ereignisse herzuleiten, den Test durchzuführen und sein Ergebnis zu interpretieren.
das Cox-Proportional-Hazards-Modell zu erklären, in R zu implementieren und Hazard Ratios fachgerecht zu interpretieren.
die zentralen Modellannahmen des Cox-Modells mithilfe von Schoenfeld- und Martingale-Residuen zu prüfen und diagnostische Auffälligkeiten zu beurteilen.
die Modellgüte anhand der Konkordanz zu beurteilen und aus den Ergebnissen begründete Handlungsempfehlungen zur Kundenabwanderung abzuleiten.

Theoretische Grundlagen der Survival Analysis

Die Survival Analysis, auch als Überlebenszeitanalyse bekannt, beschäftigt sich mit der Analyse der Zeit bis zum Eintreten eines bestimmten Ereignisses. Im Kontext der Kundenabwanderung betrachten wir die “Überlebenszeit” eines Kunden, also die Zeitspanne von der Akquisition bis zur Kündigung eines Vertrags.

Die Überlebensfunktion (Survival Function)

Die Überlebensfunktion gibt die Wahrscheinlichkeit an, dass ein Kunde länger als den Zeitpunkt $t$ “überlebt”:

\[S(t) = P(T > t)\]

wobei $T$ die Überlebenszeit (Zufallsvariable) ist. Die Survival Function gibt also die Wahrscheinlichkeit an, dass das interessierende Ereignis (z. B. Kündigung, Tod, …) nicht bis zum Zeitpunkt $t$ eintritt.

In der Praxis ist jedoch die Abschätzung dieser Wahrscheinlichkeit besonders durch sogenannte zensierte Daten schwierig, d. h. wenn beispielsweise ein Arbeitnehmer das Unternehmen aufgrund von Krankheit, Todesfall oder Wegzug ins Ausland verlässt – also aus einem Grund, der in der Studie nicht betrachtet wird. Bei medizinischen Studien wäre dies etwa der Ausstieg einer Patientin aus der Studie. Wie können wir trotz dieser “Zensierung” mancher Teilnehmender zu zuverlässigen Schätzungen gelangen?

Zensieren

In vielen praktischen Situationen (z. B. klinische Studien, Laufzeitmessungen, Kunden”lebenszeiten”) liegen zensierte Daten vor. Das heißt, für manche Probanden, Kunden oder auch Bauteile hat man nur den Zeitpunkt, bis zu dem sie beobachtet wurden (und innerhalb dieser Beobachtungszeit ist das Ereignis nicht eingetreten). Man kennt also ihre exakte „Lebensdauer“ nicht – nur, dass sie wenigstens bis zum Ende der Beobachtung durchgehalten haben. Wir unterscheiden:

Rechtszensierung (right-censoring): Wenn die Studie endet oder die Person aus anderen Gründen nicht mehr beobachtet wird, ohne dass das Ereignis aufgetreten ist.
Linkszensierung (left-censoring): Seltenere Variante, bei der man nur weiß, dass das Ereignis bereits vor der ersten Beobachtung eingetreten sein muss. (In der Praxis weniger typisch für Kaplan-Meier, aber theoretisch möglich.)
Intervallzensierung (interval-censoring): Bei Intervallzensierung weiß man nur, dass das Ereignis in einem gewissen Zeitintervall passiert sein muss, aber nicht genau, wann. Ein praktisches Beispiel wäre, dass ein Sensor die Funktion eines Bauteils zu bestimmten Zeitpunkten misst. Wenn ein Bauteil ausfällt, kann man in der Regel nur sagen, dass es zwischen zwei Messungen ausgefallen ist.

Der Kaplan-Meier-Schätzer schafft es nun, diese teilweise fehlende Information geeignet zu verarbeiten, indem er alle Personen (bzw. Objekte) bis zum Zeitpunkt, an dem wir sie verlieren, berücksichtigt und danach herausnimmt. Dies sehen wir uns jetzt an.

Kaplan-Meier-Schätzer

Statt eine globale Schätzung für die gesamte Kurve auf einmal durchzuführen, baut der Kaplan-Meier-Schätzer die Überlebenswahrscheinlichkeit schrittweise auf. Man stellt sich vor, dass zu Beginn alle Teilnehmenden (oder Systeme) noch „am Leben“ sind. Dann verfolgt man, was zu den Zeitpunkten der Ereignisse passiert.

Kernidee: Kette bedingter Wahrscheinlichkeiten

Der Kaplan-Meier-Schätzer entsteht aus einer Kette bedingter Überlebenswahrscheinlichkeiten.
Für Ereigniszeitpunkte $t_1 < t_2 < \cdots < t_k$ gilt konzeptionell:

\[P(T > t_k) = P(T > t_1)\, P(T > t_2 \mid T > t_1)\, \cdots \, P(T > t_k \mid T > t_{k-1}).\]

Das heißt: Wir „überleben“ bis $t_k$, indem wir (i) bis $t_1$ überleben und (ii) danach unter der Bedingung, dass wir bis zum vorherigen Zeitpunkt überlebt haben, auch den nächsten Zeitpunkt überleben, usw.

Sortieren
Man sortiert alle Ereigniszeiten in aufsteigender Reihenfolge: \[t_1 < t_2 < \cdots < t_k.\]
Zählen
Für jeden Ereigniszeitpunkt $t_i$ bestimmt man:
- $d_i$: Anzahl der Ereignisse, die genau zum Zeitpunkt $t_i$ auftreten.
- $n_i$: Anzahl der Personen (bzw. Geräte), die kurz vor $t_i$ noch beobachtet werden (also noch nicht ausgefallen und nicht zensiert sind). Dieses Set nennt man das Risiko-Set.
Multiplikative Aktualisierung
Der Kaplan-Meier-Schätzer ist die Produktform: \[\hat{S}(t) \;=\; \prod_{i : t_i \leq t} \Bigl(1 - \frac{d_i}{n_i}\Bigr).\]
- Idee: $\bigl(1 - d_i/n_i\bigr)$ ist die geschätzte bedingte Wahrscheinlichkeit, den Zeitpunkt $t_i$ zu „überleben“, gegeben, dass man kurz davor noch „at risk“ war.
- Durch die Multiplikation all dieser Faktoren für alle Zeitpunkte $t_i$ bis $t$ erhält man die gesamte Überlebenswahrscheinlichkeit bis zu $t$.
Schrittfunktion
Graphisch ergibt sich eine Treppenfunktion, die nur dann einen Sprung nach unten macht, wenn tatsächlich Ereignisse auftreten. Zwischen den Ereignissen bleibt sie konstant.

Wir hatten die zensierten Daten als besondere Herausforderung definiert. Wie genau wird die Zensur nun beim Kaplan-Meier-Schätzer berücksichtigt?

Im Schätzer werden zensierte Datenwerte durch $n_i$ berücksichtigt, d. h. in $\hat{S}(t)$ finden nur solche Kunden, Probanden, Bauteile etc. Berücksichtigung, die zum Zeitpunkt $t_i$ “at risk” waren. Damit fallen alle heraus, die z. B. aus der Studie aussteigen.

\[\hat{S}(t) \;=\; \prod_{i : t_i \leq t} \Bigl(1 - \frac{d_i}{n_i}\Bigr).\] Graphisch sehen wir dann Überlebenskurven, wie beispielsweise diese:

Alle Überlebenskurven fallen ab, aber wir können mit dem Kaplan-Meier-Schätzer Unterschiede zwischen Gruppen sehen. Genau für diese Unterschiede interessieren wir uns bei der Überlebenszeitanalyse. Ob jedoch die Unterschiede statistisch signifikant sind, müssen wir im nächsten Schritt untersuchen.

Der Log-Rank-Test

Die Kaplan-Meier-Kurven beschreiben das Überlebensverhalten der Gruppen. Sichtbare Abstände zwischen den Kurven können jedoch auch durch Zufall entstehen. Der Log-Rank-Test beantwortet deshalb die Frage:

Sind die über den gesamten Beobachtungszeitraum auftretenden Unterschiede grösser, als wir sie bei gleichem Ereignisrisiko in allen Gruppen erwarten würden?

Nullhypothese

Nullhypothese des Log-Rank-Tests

Unter der Nullhypothese besitzen alle Gruppen zu jedem Zeitpunkt dieselbe Hazardfunktion:

\[H_0:\quad h_1(t)=h_2(t)=\dots=h_g(t) \qquad \text{für alle }t\geq 0.\]

Bei identischem Ausgangswert $S_j(0)=1$ ist dies gleichbedeutend mit identischen Überlebensfunktionen:

\[S_1(t)=S_2(t)=\dots=S_g(t).\]

Der Test betrachtet ausschliesslich Zeitpunkte, an denen mindestens ein Ereignis beobachtet wird. Jeder Ereigniszeitpunkt ist dabei ein kleines bedingtes Zufallsexperiment: Aus welcher Gruppe stammen die Ereignisse, gegeben das aktuelle Risiko-Set?

Das Zufallsmodell an einem Ereigniszeitpunkt

Wir betrachten zunächst zwei Gruppen und einen Ereigniszeitpunkt $d_k$. Unmittelbar vor $d_k$ seien

$R_k$ das Risiko-Set,
$r_k=|R_k|$ die Gesamtzahl der Personen im Risiko-Set,
$r_{1k}$ und $r_{2k}$ die entsprechenden Anzahlen in Gruppe 1 und 2,
$q_k$ die Gesamtzahl der Ereignisse bei $d_k$,
$O_{1k}$ die tatsächlich beobachtete Anzahl dieser Ereignisse in Gruppe 1.

$O_{1k}$ ist die Zufallsvariable, deren Erwartungswert wir benötigen.

Unter $H_0$ besitzen alle Personen im Risiko-Set dieselbe momentane Ereignisrate. Bedingt auf $R_k$ und darauf, dass bei $d_k$ insgesamt $q_k$ Ereignisse auftreten, sind die $q_k$ betroffenen Personen daher wie eine Auswahl ohne Zurücklegen aus den $r_k$ Personen des Risiko-Sets. Folglich gilt

\[O_{1k}\mid R_k,q_k \sim \operatorname{Hypergeom}(r_k,r_{1k},q_k).\]

Herleitung des Erwartungswerts

Für jede Person $i$ aus Gruppe 1, die sich in $R_k$ befindet, definieren wir den Indikator

\[I_{ik}= \begin{cases} 1, & \text{wenn Person }i\text{ zu den }q_k\text{ Ereignissen gehört},\\ 0, & \text{sonst}. \end{cases}\]

Dann ist

\[O_{1k}=\sum_{i\in R_k\cap G_1} I_{ik}.\]

Unter $H_0$ hat jede der $r_k$ Personen dieselbe bedingte Wahrscheinlichkeit, zu den $q_k$ Ereignissen zu gehören:

\[P(I_{ik}=1\mid R_k,q_k)=\frac{q_k}{r_k}.\]

Mit der Linearität des Erwartungswerts folgt

\[\begin{aligned} E(O_{1k}\mid R_k,q_k) &=E\!\left(\sum_{i\in R_k\cap G_1}I_{ik}\,\middle|\,R_k,q_k\right)\\ &=\sum_{i\in R_k\cap G_1}E(I_{ik}\mid R_k,q_k)\\ &=\sum_{i\in R_k\cap G_1}\frac{q_k}{r_k}\\ &=r_{1k}\frac{q_k}{r_k}. \end{aligned}\]

Damit erhalten wir

\[\boxed{E_{1k}=\frac{r_{1k}}{r_k}\,q_k.}\]

Der oft verwendete Satz „Die Ereignisse verteilen sich proportional zur Zusammensetzung des Risiko-Sets“ ist somit keine zusätzliche Regel, sondern die Interpretation dieses bedingten Erwartungswerts.

Varianz am Ereigniszeitpunkt

Da $O_{1k}$ hypergeometrisch verteilt ist, gilt

\[\begin{aligned} \operatorname{Var}(O_{1k}\mid R_k,q_k) &=q_k\frac{r_{1k}}{r_k} \left(1-\frac{r_{1k}}{r_k}\right) \frac{r_k-q_k}{r_k-1}\\[2mm] &=\frac{r_{1k}r_{2k}q_k(r_k-q_k)} {r_k^2(r_k-1)}. \end{aligned}\]

Der letzte Faktor ist die Endlichkeitskorrektur der Ziehung ohne Zurücklegen. Bei nur einem Ereignis ($q_k=1$) reduziert sich die Situation auf die Frage, ob dieses eine Ereignis aus Gruppe 1 oder aus Gruppe 2 stammt.

Aggregation über alle Ereigniszeitpunkte

An jedem Ereigniszeitpunkt erhalten wir eine Abweichung

\[O_{1k}-E_{1k}.\]

Positive und negative Abweichungen werden zunächst mit ihrem Vorzeichen über die Zeit aufsummiert:

\[U=\sum_{k=1}^{K}(O_{1k}-E_{1k}).\]

Unter $H_0$ erwarten wir $U\approx 0$. Die zugehörige Varianz ist

\[V=\sum_{k=1}^{K}\operatorname{Var}(O_{1k}\mid R_k,q_k).\]

Die standardisierte Log-Rank-Teststatistik für zwei Gruppen lautet deshalb

\[Z=\frac{U}{\sqrt V},\]

beziehungsweise in der üblichen Chi-Quadrat-Form

\[\boxed{\chi^2=\frac{U^2}{V}\sim\chi^2_1}\]

asymptotisch unter $H_0$.

Warum nicht jede Abweichung einzeln quadrieren?

Der Log-Rank-Test sucht nach einer über die Zeit konsistenten Richtung der Abweichungen. Würden wir $(O_{1k}-E_{1k})^2$ für jeden Zeitpunkt separat summieren, gingen die Vorzeichen verloren. Auch zufällige Wechsel zwischen „zu vielen“ und „zu wenigen“ Ereignissen würden sich dann künstlich verstärken. Beim Log-Rank-Test wird deshalb zuerst die vorzeichenbehaftete Gesamtabweichung $U$ gebildet und erst danach quadriert.

Ein kleiner p-Wert spricht gegen die Nullhypothese gleicher Hazard- beziehungsweise Überlebensfunktionen. Der Test zeigt jedoch nicht, wann sich die Gruppen unterscheiden oder wie gross der Unterschied zu einem bestimmten Zeitpunkt ist. Dafür betrachten wir zusätzlich die Kaplan-Meier-Kurven und geeignete Effektmasse.

Für Interessierte: mehr als zwei Gruppen

Für $g$ Gruppen bilden wir einen Vektor der aufsummierten Abweichungen

\[\mathbf U=(U_1,\dots,U_{g-1})^\top\]

und die zugehörige Kovarianzmatrix $\mathbf V$. Eine Gruppe ist redundant, weil sich die Abweichungen über alle Gruppen zu null summieren. Die Teststatistik lautet

\[\chi^2=\mathbf U^\top\mathbf V^{-1}\mathbf U \sim\chi^2_{g-1}.\]

Damit wird zugleich verständlich, weshalb der Test bei $g$ Gruppen $g-1$ Freiheitsgrade besitzt.

Schliessen wir das Theoriekapitel mit zwei Verständnisfragen:

Verständnisfrage 1: Kaplan–Meier-Kurve – was wird eigentlich geschätzt?

Eine Kaplan–Meier-Kurve fällt nur an Zeitpunkten, an denen Ereignisse auftreten. Zensierte Beobachtungen führen zwar zu einer kleineren Risikogruppe, aber nicht zu einem „Sprung“ der Kurve.

Frage: Warum ist das so? Formuliere die Idee in 2–3 Sätzen (Stichworte: Risk-Set, bedingte Überlebenswahrscheinlichkeit, Zensierung).

Antwort anzeigen

Die Kaplan–Meier-Schätzung basiert auf dem Produkt von bedingten Überlebenswahrscheinlichkeiten: Zu jedem Ereigniszeitpunkt $d_k$ wird geschätzt, welcher Anteil der Personen im Risk-Set das Ereignis nicht erleidet.
Zensierung bedeutet: Wir wissen nur, dass die Person bis zu diesem Zeitpunkt überlebt hat, danach aber nicht mehr — daher wird sie ab diesem Zeitpunkt aus dem Risk-Set entfernt. Es passiert aber kein Ereignis, also gibt es keinen Grund für einen Sprung in der Überlebensfunktion. Sprünge entstehen nur durch beobachtete Ereignisse.

Verständnisfrage 2: Log-Rank-Test – Erwartung wovon?

Am Ereigniszeitpunkt $d_k$ wird

\[E_{1k}=\frac{r_{1k}}{r_k}\,q_k\]

berechnet.

Frage: Welche Zufallsvariable besitzt diesen Erwartungswert, auf welche Informationen wird bedingt und welche Annahme der Nullhypothese begründet die zugrunde liegende Verteilung?

Antwort anzeigen

Die Zufallsvariable $O_{1k}$ zählt, wie viele der insgesamt $q_k$ Ereignisse bei $d_k$ aus Gruppe 1 stammen. Wir bedingen auf das aktuelle Risiko-Set $R_k$ und auf die Gesamtzahl $q_k$ der Ereignisse. Unter $H_0$ haben alle Personen im Risiko-Set dieselbe Hazard. Deshalb entspricht die Gruppenzuordnung der $q_k$ Ereignisse einer Auswahl ohne Zurücklegen, also einer hypergeometrischen Verteilung. Deren Erwartungswert ist

\[E(O_{1k}\mid R_k,q_k)=q_k\frac{r_{1k}}{r_k}.\]

Damit sind wir nun bereit, uns ein konkretes Beispiel anzusehen.

Praktisches Beispiel implementieren

Wir simulieren Daten für die Abonnementlaufzeit von Kunden eines Telekommunikationsanbieters. Die Ereigniszeiten werden dabei direkt aus einem proportionalen Hazard-Modell erzeugt. Dadurch sind Ereigniszeit, Kovariaten und Zensierung statistisch konsistent miteinander verknüpft. Die Kunden unterscheiden sich im Vertragstyp, im Serviceplan und in den monatlichen Kosten.

Code

set.seed(123)
n <- 1000

# Kundenmerkmale
data <- data.frame(
  customer_id = 1:n,
  contract_type = sample(
    c("Monthly", "Annual"), n, replace = TRUE,
    prob = c(0.7, 0.3)
  ),
  service_plan = sample(
    c("Basic", "Standard", "Premium"), n, replace = TRUE
  ),
  monthly_charges = round(runif(n, 30, 150), 2)
)

# Linearer Prädiktor des proportionalen Hazard-Modells
charges_z <- as.numeric(scale(data$monthly_charges))
linear_predictor <- with(data,
  log(1.7)  * (contract_type == "Monthly") +
  log(0.95) * (service_plan == "Standard") +
  log(0.65) * (service_plan == "Premium") +
  log(1.25) * charges_z
)

# Ereigniszeiten aus einem PH-Modell mit konstanter Baseline-Hazard
baseline_hazard <- 0.035
event_time <- rexp(
  n,
  rate = baseline_hazard * exp(linear_predictor)
)

# Unabhängige Zensierung und administratives Studienende nach 36 Monaten
censor_time <- runif(n, min = 12, max = 36)
administrative_end <- 36

observed_time <- pmin(event_time, censor_time, administrative_end)
data$churn <- as.integer(
  event_time <= censor_time & event_time <= administrative_end
)
data$tenure <- pmax(round(observed_time, 1), 0.1)

Kaplan-Meier-Analyse

Code

# Survival-Objekt
surv_obj <- Surv(time = data$tenure, event = data$churn)

# Gesamte Überlebenskurve
km_fit <- survfit(surv_obj ~ 1)

# Stratifiziert nach Vertragstyp
km_contract <- survfit(surv_obj ~ contract_type, data = data)

# Plots
ggsurvplot(
  km_fit,
  data = data,
  title = "Gesamte Überlebenskurve",
  xlab = "Zeit (Monate)",
  ylab = "Überlebenswahrscheinlichkeit",
  conf.int = TRUE,
  risk.table = TRUE
)

Code

ggsurvplot(
  km_contract,
  data = data,
  title = "Überlebenskurven nach Vertragstyp",
  xlab = "Zeit (Monate)",
  ylab = "Überlebenswahrscheinlichkeit",
  conf.int = TRUE,
  pval = TRUE
)

Die erste Überlebenszeitkurve (ohne Strata) umfasst alle Kunden des Unternehmens. Beim zweiten Plot wurden zwei Gruppen gebildet, je nach monatlicher oder jährlicher Laufzeit. Wir beobachten bis ca. 12 Monate kaum Unterschiede, jedoch erscheint danach die Anzahl Kündigungen bei Kunden mit monatlicher Laufzeit (und damit mehr Flexibilität bei der Kündigung) höher auszufallen als bei denen mit jährlicher Frist. Ob die Unterschiede jedoch statistisch signifikant sind, muss jetzt der Log-Rank-Test beantworten.

Log-Rank-Test

Code

# Log-Rank-Test für den Vertragstyp
log_rank_test <- survdiff(surv_obj ~ contract_type, data = data)
print(log_rank_test)

Call:
survdiff(formula = surv_obj ~ contract_type, data = data)

                        N Observed Expected (O-E)^2/E (O-E)^2/V
contract_type=Annual  295      157      214      15.1        23
contract_type=Monthly 705      477      420       7.7        23

 Chisq= 23  on 1 degrees of freedom, p= 2e-06

Code

# survdiff() gibt die Teststatistik, aber nicht direkt den p-Wert aus
log_rank_df <- length(log_rank_test$n) - 1
p_logrank <- pchisq(log_rank_test$chisq, df = log_rank_df, lower.tail = FALSE)
p_logrank

[1] 1.619061e-06

Der berechnete p-Wert beträgt 1.62e-06. Wir vergleichen ihn mit dem vorab festgelegten Signifikanzniveau $\alpha$. Ist $p<\alpha$, verwerfen wir die Nullhypothese gleicher Hazardfunktionen. Die Kaplan-Meier-Kurven benötigen wir trotzdem weiterhin: Der Test verdichtet den gesamten Verlauf auf eine Testentscheidung, während die Kurven zeigen, in welcher Richtung und in welchen Zeitbereichen die Unterschiede auftreten.

Schauen wir uns die Kaplan-Meier-Plots nochmals genauer an. Neben den Unterschieden der abfallenden Treppenfunktionen sind hier noch die Bänder der Konfidenzintervalle farbig eingezeichnet. Diese werden mit zunehmender Zeit grösser, da die Unsicherheit der Schätzung zunimmt. Dies vor allem, weil die Anzahl der Teilnehmenden immer kleiner wird, d. h. wir beziehen Informationen aus einer immer kleiner werdenden Gruppe. Dies tritt ggf. auch schon früher auf, wenn wir z. B. eine kleine Gruppe mit einer großen Gruppe vergleichen.

Exkurs: Bootstrap bei zensierten Daten

Bei kleinen Risiko-Sets werden Kaplan–Meier-Schätzungen und asymptotische Konfidenzintervalle zunehmend unsicher. Bootstrapverfahren können bei rechtszensierten Daten zur Schätzung von Standardfehlern oder zur Konstruktion von Konfidenzbändern eingesetzt werden. Sie ersetzen jedoch keine ausreichende Zahl beobachteter Ereignisse und sind deshalb nicht automatisch eine Lösung für sehr kleine Stichproben. Zudem muss das Resampling die Zensierungsstruktur angemessen berücksichtigen (Efron, 1981; Akritas, 1986).

Dieser Ansatz geht über die Zielsetzung des Tutorials hinaus.

Fazit: Was lernen wir?

Mit dem Log-Rank-Test und den Kaplan-Meier-Schätzern können wir paarweise die Überlebenszeitkurven gut vergleichen und Unterschiede evaluieren. Jedoch ist das eher explorativer Natur. Wenn wir mehrere Variablen haben, können wir auf diese Art kein komplettes Modell zusammenstellen und evaluieren, wie wir es z. B. von der multiplen linearen Regression oder der logistischen Regression gewohnt sind.

Um dies zu können, wurde das Cox-Regressionsmodell (oder auch Cox Proportional Hazards Model) entwickelt, welches wir uns im nächsten Abschnitt ansehen.

Cox-Regressionsmodell bzw. Cox Proportional Hazards Model

In der Überlebenszeitanalyse interessieren wir uns für den Einfluss unserer Variablen oder auch Kovariaten auf die Überlebenszeit eines Objekts. Beim Kündigungsverhalten unserer Abokunden fragen wir uns, wie groß der Einfluss (“Hazard”) der erhobenen Variablen, wie die Wahl des Abotyps oder die Höhe der monatlichen Kosten ist. Wie bei der multiplen linearen Regression möchten wir wissen, ob wir die einflussreichsten Variablen/Kovariaten identifizieren können, um ein möglichst einfaches Modell zu entwickeln.

Das Ziel der Überlebenszeitanalyse ist nicht nur eine oder mehrere Kurve(n) zu zeichnen, sondern eine Analysefrage zu beantworten:

Welche Merkmale erhöhen oder senken die Kündigungsrate — und um wie viel?

Dafür brauchen wir eine Grösse, die drei Dinge gleichzeitig kann:

Sie funktioniert auch mit Zensierung (Kunden, die am Ende der Beobachtung noch nicht gekündigt haben).
Sie ist eine Momentaufnahme: Wir wollen das Risiko jetzt vergleichen (nicht nur „irgendwann“).
Sie erlaubt einen fairen Vergleich zwischen Kunden, die unterschiedlich lange schon dabei sind (tenure).

Einige werden jetzt einwenden, dass die Momentaufnahme mich doch gar nicht interessiert um am Ende eine Businessentscheidung zu treffen und z.B. das Abo-Modell ändere. Das ist natürlich völlig korrekt, wenn ich aber einen grösseren Zeitraum betrachte -“wie hoch ist das Risiko, dass ein bestimmter Kunde innerhalb der nächsten 12 Monate kündigt”- vermische ich den Einfluss mehrerer Effekte (z.B. Saisonalität, Länge der Kundenbeziehung). In diesem Zeitraum habe ich alle möglichen Kunden die kündigen. Kann ich jemanden, der nach 4 Monaten kündigt mit jemanden der nach 10 Monaten kündigt, vergleichen? Wahrscheinlich nicht. Betrachte ich es jedoch zeitlokal vergleiche ich nicht Äpfel mit Birnen, sondern Kunden, die alle zu diesem Zeitpunkt at risk waren.

Genau das leistet die Hazard oder Hazard-Funktion.

Was misst die Hazard?

Anschaulich: „Unter den Kunden, die bis heute noch nicht gekündigt haben: Wie gross ist das momentane Kündigungsrisiko in der nächsten kurzen Zeitspanne?“

Damit wird das Risiko vergleichbar: Wir vergleichen Kunden unter derselben Bedingung („noch im Bestand“).

Die Hazard-Funktion

Die Hazard-Funktion beschreibt das momentane Risiko einer Kündigung unter der Bedingung $T > t$:

\[\begin{aligned} h(t) &= \lim_{\Delta t \to 0} \frac{P(t < T \le t + \Delta t \mid T > t)}{\Delta t} \\ &= \lim_{\Delta t \to 0} \frac{P(t < T \le t + \Delta t \ \cap\ T > t)}{\Delta t \cdot P(T>t)} \\ &= \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t)}{\Delta t }\cdot \frac{1}{P(T>t)} \quad (\text{da aus } t \le T < t+\Delta t \text{ automatisch } T>t \text{ folgt})\\ &= f(t) \cdot \frac{1}{S(t)} \quad (\text{da } S(t)=P(T>t)). \end{aligned}\]

Daraus folgt dann, dass wir $f(t)$ aus $h(t)$ und $S(t)$ berechnen können: $f(t)=h(t)\cdot S(t)$. In der Praxis kann man zwar öfters auch die vereinfachende Annahme machen, dass $h(t) = \text{const.}$ gilt, wir wollen aber gleich zum flexibleren proportional hazard model weitergehen.

Das Cox Proportional Hazards Model

Die zentrale Vergleichsgröße: Hazard Ratio

In der Praxis wollen wir selten das absolute Kündigungsrisiko zu einem Zeitpunkt $t$ wissen. Viel wichtiger ist der Vergleich:

Hat Person 1 zum Zeitpunkt $t$ ein höheres oder niedrigeres Kündigungsrisiko als Person 2 – und um wie viel?

Dazu betrachten wir die Hazard Ratio (HR) zweier Personen mit Kovariaten $x^{(1)}$ und $x^{(2)}$:

\[\mathrm{HR}(t) = \frac{h\!\left(t\mid x^{(1)}\right)}{h\!\left(t\mid x^{(2)}\right)}.\]

Interpretation:

$\mathrm{HR}(t)=1$: gleiches momentanes Risiko,
$\mathrm{HR}(t)=2$: Person 1 hat das doppelte momentane Risiko wie Person 2,
$\mathrm{HR}(t)=0.5$: Person 1 hat das halbe momentane Risiko.

Der Vorteil: Die HR ist eine relative Grösse. Sie bleibt auch dann interpretierbar, wenn das Basisrisiko über die Zeit stark schwankt (Saisonalität (“alle Arten von Fitnessstudio-Effekt im Januar”), Servicegebühren, die nur zwei Mal im Jahr in Rechnung gestellt werden etc. ).

Ein Modell, das das Arbeiten mit Hazard Ratios einfach macht: Der Separationsansatz

Damit wir den Einfluss von Kovariaten systematisch beschreiben können, wählen wir ein Modell, das die Hazard in einen Zeitanteil und einen Kovariatenanteil trennt:

\[h(t\mid x)=h_0(t)\,g(x).\]

$h_0(t)$ ist das zeitabhängige Basisrisiko (wie hoch ist das Kündigungsrisiko grundsätzlich über die Zeit?),
$g(x)$ beschreibt den relativen Einfluss der Kovariaten.

Der entscheidende Punkt: Für zwei Personen kürzt sich das Basisrisiko in der Hazard Ratio heraus:

\[\frac{h(t\mid x^{(1)})}{h(t\mid x^{(2)})} = \frac{h_0(t)g(x^{(1)})}{h_0(t)g(x^{(2)})} = \frac{g(x^{(1)})}{g(x^{(2)})}\]

Damit hängt die Hazard Ratio nicht mehr von $t$ ab — genau das ist die Proportional-Hazards-Idee.

Cox Proportional Hazards: Wahl von $g(x)$

Wir wählen eine stets positive und gut interpretierbare Form:

\[g(x)=\exp(\boldsymbol{\beta^\top x}),\] Anm.: $\boldsymbol{\beta^T} \boldsymbol{x}$ ist ein sog. Skalarprodukt aus dem Koeffizientenvektor $\boldsymbol{\beta}$ und dem Vektor der Kovariaten $\boldsymbol{x}$ und die kompakte Schreibweise von:

\[\boldsymbol{\beta}^{\top}\boldsymbol{x} = \begin{pmatrix} \beta_1 & \beta_2 & \cdots & \beta_p \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_p \end{pmatrix} = \beta_1x_1+\beta_2x_2+\cdots+\beta_px_p.\]

und erhalten das Cox-Modell

\[h(t|x)=h_0(t)\exp(\beta^\top x).\]

Dann gilt für die Hazard Ratio:

\[\frac{h(t|x^{(1)})}{h(t|x^{(2)})} = \exp\!\left(\beta^\top(x^{(1)}-x^{(2)})\right),\]

und insbesondere für eine Erhöhung von $x_j$ um 1 (bei konstanten übrigen Variablen):

\[\mathrm{HR}=\exp(\beta_j).\]

Cox-Modell und Proportional Hazards – anschaulich am Kundenbeispiel

In unseren Daten beschreiben die Kovariaten $\boldsymbol{x}$ Eigenschaften eines Kunden:

contract_type (Monthly vs. Annual)
service_plan (Basic / Standard / Premium)
monthly_charges (Monatskosten)

Die zentrale Frage lautet: Welche dieser Kovariaten erhöhen oder senken die Kündigungsrate?

Was bedeutet Proportional Hazards konkret?

Stellen wir uns zwei Kunden vor, die sich nur in einer Kovariaten unterscheiden:

Kunde A: contract_type = Monthly
Kunde B: contract_type = Annual

Wenn „Monthly“ grundsätzlich riskanter ist, dann erwarten wir:

\[\mathrm{HR}(t)=\frac{h(t|\text{Monthly})}{h(t|\text{Annual})} > 1.\]

Proportional Hazards bedeutet jetzt etwas sehr Spezifisches:

Dieses Verhältnis (die Hazard Ratio) ist über die Zeit konstant.

Also: Der Monthly Kunde besitzt z. B. zu jedem Zeitpunkt ein 1.8-mal so hohes Risiko wie der Annual-Kunde, unabhängig davon, ob wir am Anfang der Kundenbeziehung sind oder nach vielen Monaten.

Wichtig

Das Basisrisiko darf sich trotzdem über die Zeit ändern (z. B. viele Kündigungen zu Beginn, später weniger). Proportional heißt nur: die Kurven sind vertikal skaliert, nicht identisch!

Wie kommen wir an die Koeffizienten $\beta$?

Für Kunde $i$ beobachten wir

$\tilde t_i$: Ereignis- oder Zensierungszeit,
$\delta_i\in\{0,1\}$: Ereignisindikator,
$x_i$: Kovariatenvektor.

Ist $\delta_i=1$, kennen wir den Ereigniszeitpunkt. Ist $\delta_i=0$, wissen wir nur, dass die tatsächliche Ereigniszeit grösser als $\tilde t_i$ ist.

Vollständige Likelihood

Der Likelihood-Beitrag einer Beobachtung lautet

\[L_i=f(\tilde t_i\mid x_i)^{\delta_i} S(\tilde t_i\mid x_i)^{1-\delta_i}.\]

Mit $f(t\mid x)=h(t\mid x)S(t\mid x)$ erhalten wir äquivalent

\[L_i=h(\tilde t_i\mid x_i)^{\delta_i} S(\tilde t_i\mid x_i).\]

Im Cox-Modell enthält diese vollständige Likelihood neben $\beta$ auch die unbekannte Funktion $h_0(t)$. Cox’ entscheidende Idee besteht deshalb darin, nur die Information zu verwenden, die für den relativen Kovariateneffekt relevant ist.

Bedingte Wahrscheinlichkeit innerhalb des Risiko-Sets

Betrachten wir einen Ereigniszeitpunkt $d_k$, an dem genau eine Person $i(k)$ kündigt. Das Risiko-Set $R(d_k)$ enthält alle Personen, die unmittelbar davor noch beobachtet und ereignisfrei sind.

Bedingt darauf, dass bei $d_k$ ein Ereignis eintritt, ist die Wahrscheinlichkeit, dass gerade Person $i(k)$ betroffen ist,

\[P\!\left(i(k)\text{ kündigt}\mid R(d_k),\text{ ein Ereignis bei }d_k\right) = \frac{h(d_k\mid x_{i(k)})} {\sum_{j\in R(d_k)}h(d_k\mid x_j)}.\]

Setzen wir

\[h(t\mid x)=h_0(t)\exp(\beta^\top x)\]

ein, kürzt sich das für alle Personen gleiche Basisrisiko heraus:

\[\frac{h_0(d_k)\exp(\beta^\top x_{i(k)})} {\sum_{j\in R(d_k)}h_0(d_k)\exp(\beta^\top x_j)} = \frac{\exp(\beta^\top x_{i(k)})} {\sum_{j\in R(d_k)}\exp(\beta^\top x_j)}.\]

Das Produkt über alle Ereigniszeitpunkte ergibt die partielle Likelihood:

\[L_p(\beta)= \prod_{k=1}^{K} \frac{\exp(\beta^\top x_{i(k)})} {\sum_{j\in R(d_k)}\exp(\beta^\top x_j)}.\]

Wir maximieren ihre Logarithmusfunktion numerisch:

\[\ell_p(\beta) = \sum_{k=1}^{K} \left[ \beta^\top x_{i(k)} - \log\!\left( \sum_{j\in R(d_k)}\exp(\beta^\top x_j) \right) \right].\]

Was ist an der Likelihood partiell?

Sie ist keine vollständige Wahrscheinlichkeitsbeschreibung der beobachteten Zeiten. Sie verwendet nur die bedingte Information darüber, wer innerhalb des jeweiligen Risiko-Sets das Ereignis erlebt. Genau dadurch kann $\beta$ geschätzt werden, ohne $h_0(t)$ parametrisch festzulegen.

Bei mehreren Ereignissen zum selben Zeitpunkt verwendet Software Näherungen beziehungsweise Tie-Verfahren wie Efron oder Breslow.

Praktische Implementierung unseres Beispiels

Cox Proportional Hazards Model

Code

# Modell-Fit
cox_model <- coxph(
  surv_obj ~ contract_type + service_plan + monthly_charges,
  data = data
)

# Zusammenfassung
print(summary(cox_model))

Call:
coxph(formula = surv_obj ~ contract_type + service_plan + monthly_charges, 
    data = data)

  n= 1000, number of events= 634 

                          coef exp(coef)  se(coef)      z Pr(>|z|)    
contract_typeMonthly  0.490424  1.633009  0.092789  5.285 1.25e-07 ***
service_planPremium  -0.504699  0.603687  0.100008 -5.047 4.50e-07 ***
service_planStandard -0.113130  0.893035  0.094281 -1.200     0.23    
monthly_charges       0.008282  1.008316  0.001170  7.076 1.48e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                     exp(coef) exp(-coef) lower .95 upper .95
contract_typeMonthly    1.6330     0.6124    1.3615    1.9587
service_planPremium     0.6037     1.6565    0.4962    0.7344
service_planStandard    0.8930     1.1198    0.7424    1.0743
monthly_charges         1.0083     0.9918    1.0060    1.0106

Concordance= 0.618  (se = 0.012 )
Likelihood ratio test= 107.6  on 4 df,   p=<2e-16
Wald test            = 104.6  on 4 df,   p=<2e-16
Score (logrank) test = 106.2  on 4 df,   p=<2e-16

Modelldiagnostik

Wie bei anderen Regressionsmodellen müssen wir prüfen, ob die Struktur des Cox-Modells zu den Daten passt. Zwei Residuentypen beantworten dabei unterschiedliche Fragen:

Schoenfeld-Residuen: Bleibt der Effekt einer Kovariate über die Zeit konstant?
Martingale-Residuen: Ist die gewählte Funktionsform einer metrischen Kovariate angemessen?

Schoenfeld-Residuen: Ist der Effekt zeitkonstant?

Für jedes beobachtete Ereignis und jede Kovariate $j$ wird der Wert der tatsächlich kündigenden Person mit einem vom Modell erwarteten Wert im aktuellen Risiko-Set verglichen:

\[r^{(S)}_{ij} = x_{ij}-\bar x_j(t_i;\hat\beta),\]

wobei

\[\bar x_j(t_i;\hat\beta) = \frac{\sum_{\ell\in R(t_i)}x_{\ell j} \exp(\hat\beta^\top x_\ell)} {\sum_{\ell\in R(t_i)} \exp(\hat\beta^\top x_\ell)}.\]

Wenn die Proportional-Hazards-Annahme gilt, sollten die skalierten Schoenfeld-Residuen keinen systematischen Trend über die Zeit zeigen. Entscheidend ist deshalb nicht, ob die einzelnen Punkte stark streuen, sondern ob die geglättete Kurve zeitlich ansteigt, fällt oder sich systematisch krümmt.

Warum sehen die Punktwolken je nach Variablentyp anders aus?

Bei einer binären Kovariate ist der beobachtete Wert $x_{ij}$ entweder 0 oder 1. Die Residuen bilden deshalb häufig zwei Bänder. Bei einer metrischen Kovariate kann $x_{ij}$ viele Werte annehmen; die Punkte verteilen sich entsprechend über einen grösseren vertikalen Bereich. Diese unterschiedliche Geometrie ist normal. Für die PH-Annahme zählt der Zeittrend, nicht die Form der Punktwolke an sich.

Code

ph_test <- cox.zph(cox_model)
print(ph_test)

                 chisq df    p
contract_type   0.5713  1 0.45
service_plan    2.7764  2 0.25
monthly_charges 0.0128  1 0.91
GLOBAL          3.4263  4 0.49

Code

plot(ph_test)

Ein kleiner p-Wert im Test von cox.zph() beziehungsweise ein deutlicher Trend im Plot spricht gegen einen zeitkonstanten Effekt. Ein nicht signifikanter Test beweist die PH-Annahme jedoch nicht; insbesondere bei kleinen Stichproben kann die Teststärke gering sein.

Martingale-Residuen: Beobachtete minus erwartete Ereignisse

Motivation: Was ist hier überhaupt ein Residuum?

In einer linearen Regression ist die Idee unmittelbar:

\[\text{Residuum}=\text{beobachteter Wert}-\text{vorhergesagter Wert}.\]

Bei Survivaldaten gibt es jedoch keinen vollständig beobachteten Zielwert für jede Person. Für eine zensierte Person kennen wir die Ereigniszeit gerade nicht. Trotzdem liefert jede Person zwei Informationen:

ob während ihrer Beobachtungszeit ein Ereignis auftrat,
wie lange sie dem Ereignisrisiko ausgesetzt war.

Ein geeignetes Residuum muss daher die beobachtete Zahl der Ereignisse mit der bis zum Beobachtungsende vom Modell erwarteten Zahl der Ereignisse vergleichen.

Schritt 1: Der beobachtete Ereignisprozess

Für Person $i$ definieren wir den Zählprozess

\[N_i(t) = \begin{cases} 0, & \text{solange bis }t\text{ kein Ereignis beobachtet wurde},\\ 1, & \text{nach dem beobachteten Ereignis}. \end{cases}\]

Bei einer zensierten Person bleibt $N_i(t)$ während der gesamten Beobachtung gleich 0. Zusätzlich bezeichnet

\[Y_i(t)= \begin{cases} 1, & \text{wenn Person }i\text{ unmittelbar vor }t\text{ noch at risk ist},\\ 0, & \text{sonst} \end{cases}\]

ob die Person zu diesem Zeitpunkt überhaupt noch zum Risiko-Set gehört.

Schritt 2: Was erwartet das Cox-Modell in einem kleinen Zeitintervall?

Unter dem Cox-Modell besitzt Person $i$ während ihrer Risikozeit die momentane Ereignisrate

\[\lambda_i(t) = Y_i(t)h_0(t)\exp(\beta^\top x_i).\]

Für ein sehr kleines Intervall $[t,t+dt)$ bedeutet dies näherungsweise

\[E\!\left[dN_i(t)\mid\text{bisheriger Verlauf}\right] = \lambda_i(t)\,dt.\]

Die rechte Seite ist also der vom Modell erwartete kleine Zuwachs der Ereigniszahl.

Schritt 3: Erwartete Ereignisse über die gesamte Beobachtungszeit

Addieren beziehungsweise integrieren wir diese kleinen erwarteten Zuwächse bis zum Zeitpunkt $t$, erhalten wir

\[A_i(t) = \int_0^t Y_i(u)h_0(u)\exp(\beta^\top x_i)\,du.\]

$A_i(t)$ wird kumulative Intensität oder Kompensator genannt. Inhaltlich ist es die bis $t$ aufgelaufene erwartete Ereigniszahl. Für zeitlich konstante Kovariaten und den beobachteten Endzeitpunkt $\tilde t_i$ gilt

\[A_i(\tilde t_i) = H_0(\tilde t_i)\exp(\beta^\top x_i),\]

wobei $H_0(t)=\int_0^t h_0(u)\,du$ die kumulative Baseline-Hazard ist.

Schritt 4: Beobachtet minus erwartet

Nun können wir ein Residuum bilden:

\[M_i(t)=N_i(t)-A_i(t).\]

Für ein kleines zukünftiges Intervall gilt bei einem korrekt spezifizierten Modell

\[E\!\left[dM_i(t)\mid\text{bisheriger Verlauf}\right] = E[dN_i(t)\mid\text{bisheriger Verlauf}]-\lambda_i(t)dt =0.\]

Genau diese Eigenschaft gibt dem Prozess seinen Namen: $M_i(t)$ ist ein Martingal. Anschaulich besitzt die Differenz „beobachtet minus erwartet“ keinen anhand der bisherigen Information vorhersagbaren systematischen Drift. Sie kann zufällig steigen oder fallen, aber das Modell sollte keine verbleibende Richtung vorhersagen können.

Am Ende der Beobachtung ist $N_i(\tilde t_i)=\delta_i$. Nach Einsetzen der geschätzten Modellgrössen erhalten wir das Martingale-Residuum

\[\boxed{ \hat M_i = \delta_i- \hat H_0(\tilde t_i)\exp(\hat\beta^\top x_i) }\]

oder kurz

\[\hat M_i = \text{beobachtete Ereigniszahl} - \text{erwartete Ereigniszahl}.\]

Interpretation eines einzelnen Martingale-Residuums

$\hat M_i>0$: Es wurde ein Ereignis beobachtet, obwohl das Modell bis zu diesem Zeitpunkt weniger als ein Ereignis erwartet hatte. Anschaulich trat das Ereignis früher oder bei geringerem vorhergesagtem kumulativem Risiko auf als vom Modell nahegelegt.
$\hat M_i\approx 0$: Beobachtete und erwartete Ereigniszahl stimmen ungefähr überein.
$\hat M_i<0$: Das Modell hat bis zum Beobachtungsende mehr Ereignisrisiko akkumuliert, als tatsächlich beobachtet wurde. Dies tritt besonders bei zensierten Personen mit hohem vorhergesagtem Risiko oder bei später als erwartet eintretenden Ereignissen auf.

Martingale-Residuen sind nicht symmetrisch

Für eine Person kann höchstens ein Ereignis beobachtet werden. Deshalb gilt $\hat M_i\leq 1$. Nach unten sind die Residuen dagegen nicht entsprechend begrenzt. Zensierte Beobachtungen haben stets $\hat M_i=-\hat A_i(\tilde t_i)\leq0$. Die Residuen sind daher typischerweise schief verteilt und nicht wie gewöhnliche normalverteilte Regressionsfehler zu lesen.

Wie zeigen Martingale-Residuen eine falsche Funktionsform?

Nehmen wir an, der wahre Einfluss der monatlichen Kosten sei gekrümmt, wir modellierten ihn aber als linearen Term $\beta\,\texttt{monthly\_charges}$. Dann bliebe ein Teil dieses systematischen Zusammenhangs im Residuum zurück.

Eine anschauliche Diagnose besteht darin, zunächst ein Cox-Modell mit den übrigen Kovariaten, aber ohne die zu untersuchende metrische Variable zu schätzen. Anschliessend tragen wir die Martingale-Residuen gegen diese Variable auf:

Eine ungefähr geradlinige geglättete Beziehung spricht für einen linearen Term im Log-Hazard.
Eine deutliche Krümmung spricht für eine Transformation, einen polynomialen Term oder einen Spline.
Einzelne extreme Punkte können auf ungewöhnliche Beobachtungen hinweisen, sind wegen der asymmetrischen Residuenverteilung aber vorsichtig zu interpretieren.

Code

# Modell ohne monthly_charges: Die Residuen sollen die noch fehlende
# Funktionsform dieser Kovariate sichtbar machen.
cox_without_charges <- coxph(
  surv_obj ~ contract_type + service_plan,
  data = data
)

martingale_plot_data <- data.frame(
  monthly_charges = data$monthly_charges,
  martingale = residuals(cox_without_charges, type = "martingale")
)

ggplot(
  martingale_plot_data,
  aes(x = monthly_charges, y = martingale)
) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(
    title = "Martingale-Residuen und Funktionsform der Monatskosten",
    x = "Monatliche Kosten",
    y = "Martingale-Residuum"
  )

Die geglättete Kurve ist hier das zentrale Diagnoseinstrument. Sie soll nicht zwingend horizontal bei null liegen, weil monthly_charges im Diagnosemodell bewusst fehlt. Ihre Form zeigt, welche funktionale Beziehung noch in das Modell aufgenommen werden sollte. Nach Wahl einer passenden Form sollte im vollständig spezifizierten Modell keine ausgeprägte systematische Reststruktur verbleiben.

Schoenfeld und Martingale prüfen unterschiedliche Annahmen

Schoenfeld-Residuen fragen: Ändert sich der Koeffizient über die Zeit?
Martingale-Residuen fragen: Ist die Form des Kovariateneffekts richtig modelliert?

Eine Kovariate kann linear im Log-Hazard wirken, aber die PH-Annahme verletzen — oder umgekehrt einen zeitkonstanten, jedoch nichtlinearen Effekt besitzen.

Modellgüte im Cox-Modell (Analogien zu $R^2$)

In der Survival-Analyse gibt es kein einziges $R^2$ wie bei der linearen Regression. Stattdessen beurteilen wir die Güte eines Cox-Modells über Diskrimination, Signifikanz und optional Kalibrierung/Fehlermaße.

Diskrimination: Konkordanz (Harrell’s C)

Was misst C? Anteil der (vergleichbaren) Paare, in denen derjenige mit kürzerer Überlebenszeit auch höheren vorhergesagten Hazard (grösseren linearen Prädiktor) hat.
Skala: 0.5 = Zufall, 1 = perfekt. Eine Konkordanz von 0.5 entspricht einer zufälligen Rangordnung, 1 einer perfekten Rangordnung. Die praktische Bedeutung von Zwischenwerten ist stark vom Anwendungskontext und vom Grad der Zensierung abhängig.

Code

# Harrell's C aus der Modellzusammenfassung
summary(cox_model)$concordance

         C      se(C) 
0.61826379 0.01197606

Code

# Aktuelle survival-Schnittstelle
concordance(cox_model)

Call:
concordance.coxph(object = cox_model)

n= 1000 
Concordance= 0.6183 se= 0.01198
concordant discordant     tied.x     tied.y    tied.xy 
    250321     154555          7        900          0

Interpretation und praktische Empfehlungen

Zunächst bereiten wir die Koeffizienten als Hazard Ratios mit Konfidenzintervallen auf:

Code

cox_summary <- summary(cox_model)

cox_results <- data.frame(
  Variable = rownames(cox_summary$coefficients),
  Hazard_Ratio = cox_summary$conf.int[, "exp(coef)"],
  CI_unten = cox_summary$conf.int[, "lower .95"],
  CI_oben = cox_summary$conf.int[, "upper .95"],
  p_Wert = cox_summary$coefficients[, "Pr(>|z|)"],
  row.names = NULL
)

knitr::kable(
  cox_results,
  digits = 3,
  caption = "Hazard Ratios des Cox-Modells"
)

Hazard Ratios des Cox-Modells
Variable	Hazard_Ratio	CI_unten	CI_oben	p_Wert
contract_typeMonthly	1.633	1.361	1.959	0.00
service_planPremium	0.604	0.496	0.734	0.00
service_planStandard	0.893	0.742	1.074	0.23
monthly_charges	1.008	1.006	1.011	0.00

Für die Interpretation gelten folgende Regeln:

Eine Hazard Ratio grösser als 1 steht für eine höhere momentane Kündigungsrate gegenüber der Referenzkategorie beziehungsweise bei einer Erhöhung der Kovariate um eine Einheit.
Eine Hazard Ratio kleiner als 1 steht für eine niedrigere momentane Kündigungsrate.
Enthält das 95%-Konfidenzintervall den Wert 1, ist der entsprechende Effekt auf dem 5%-Niveau nicht statistisch signifikant.
Bei monthly_charges bezieht sich die Hazard Ratio auf eine Geldeinheit. Für eine fachlich besser verständliche Änderung um beispielsweise 10 Einheiten kann die Hazard Ratio als $\exp(10\hat\beta)$ berechnet werden.

Da die Daten aus einem proportionalen Hazard-Modell simuliert wurden, sollten die geschätzten Richtungen ungefähr der Datengenerierung entsprechen: Monthly besitzt gegenüber Annual eine höhere Hazard, Premium gegenüber Basic eine niedrigere Hazard und steigende Monatskosten erhöhen das Risiko. Abweichungen der geschätzten Zahlen von den vorgegebenen Simulationseffekten entstehen durch Stichprobenzufall und Zensierung.

Was dürfen wir daraus praktisch ableiten?

Das Modell identifiziert statistische Zusammenhänge, aber noch keine kausalen Wirkungen. Ein höheres Kündigungsrisiko bei monatlichen Verträgen bedeutet beispielsweise nicht automatisch, dass eine Umstellung des Vertragstyps die Kündigungen im gleichen Ausmass senken würde. Für Entscheidungen können die Ergebnisse dennoch genutzt werden, um

Kundengruppen mit erhöhtem prognostiziertem Risiko genauer zu untersuchen,
gezielte Retention-Massnahmen zu testen,
mögliche nichtlineare Preiseffekte oder zeitabhängige Vertragseffekte weiter zu modellieren,
die Prognose auf unabhängigen Daten hinsichtlich Diskrimination und Kalibrierung zu validieren.

Trennung von Inferenz und Prognose

Signifikante Koeffizienten zeigen Evidenz für Zusammenhänge im Modell. Eine brauchbare individuelle Vorhersage erfordert zusätzlich eine ausreichende Diskrimination, eine gute Kalibrierung und eine externe oder zumindest resampling-basierte Validierung.

Referenzen

Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457–481.
Mantel, N. (1966). Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemotherapy Reports, 50(3), 163–170.
Efron, B. (1981). Censored data and the bootstrap. Journal of the American Statistical Association, 76(374), 312–319.
Akritas, M. G. (1986). Bootstrapping the Kaplan–Meier estimator. Journal of the American Statistical Association, 81(396), 1032–1038.
Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society. Series B (Methodological), 34(2), 187–202.
Andersen, P. K., & Gill, R. D. (1982). Cox’s regression model for counting processes: A large sample study. The Annals of Statistics, 10(4), 1100–1120.
Schoenfeld, D. (1982). Partial residuals for the proportional hazards regression model. Biometrika, 69(1), 239–241.
Therneau, T. M., Grambsch, P. M., & Fleming, T. R. (1990). Martingale-based residuals for survival models. Biometrika, 77(1), 147–160.
Grambsch, P. M., & Therneau, T. M. (1994). Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81(3), 515–526.
Klein, J. P., & Moeschberger, M. L. (2003). Survival Analysis: Techniques for Censored and Truncated Data (2nd ed.). Springer.
Harrell, F. E., Jr. (2015). Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis (2nd ed.). Springer.

--- title: "Review-Tutorial Survival Analysis" author: - name: "Markus Geuss" affiliation: "Fernfachhochschule Schweiz" title-block-banner: ./images/ffhs-farbwelt-verlauf_01.jpg lang: de format: html: theme: cosmo toc: true toc-title: "Inhaltsverzeichnis" toc-depth: 3 toc-location: right code-fold: true code-tools: true self-contained: false css: styles.css include-in-header: text: | <style> .title, .subtitle { color: white !important; } .quarto-title-banner .quarto-title { padding-left: 160px; } .aufgabe { background-color: #eaf3fb; border: 1px solid #b8d4ea; border-radius: 6px; padding: 1rem 1.25rem; margin: 1.5rem 0; } .aufgabe summary { cursor: pointer; font-weight: 600; color: #0d6efd; padding: 0.4rem 0; } .aufgabe summary:hover { text-decoration: underline; } .aufgabe details[open] summary { margin-bottom: 0.5rem; } </style> <script> document.addEventListener('DOMContentLoaded', function () { var banner = document.querySelector('.quarto-title-banner'); if (banner) { banner.style.position = 'relative'; var logo = document.createElement('img'); logo.src = './images/FFHS_Logo.png'; logo.style.cssText = [ 'position: absolute', 'left: 0px', 'top: 50%', 'transform: translateY(-50%)', 'height: 75px', 'width: auto', 'z-index: 10', 'filter: brightness(0) invert(1)' ].join(';'); banner.insertBefore(logo, banner.firstChild); } }); </script> execute: message: false warning: false editor: visual --- ```{r} #| label: setup #| include: false knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) options(repos = c(CRAN = "https://cloud.r-project.org")) suppressPackageStartupMessages({ library(survival) library(ggplot2) library(dplyr) library(survminer) }) theme_set(theme_minimal(base_size = 12)) ``` ::: {.callout-tip title="Lernziele"} Nach dem Erarbeiten des Tutorials sind Sie in der Lage … 1. die grundlegenden Konzepte der Survival Analysis (Überlebenszeit, Überlebensfunktion, Hazard-Rate und Zensierung) zu erläutern. 2. Überlebensdaten in R aufzubereiten und mit einer Zeitvariable sowie einem Ereignisindikator korrekt abzubilden. 3. den Kaplan-Meier-Schätzer aus bedingten Überlebenswahrscheinlichkeiten zu erklären, in R anzuwenden und die geschätzten Überlebenskurven zu interpretieren. 4. die Idee des Log-Rank-Tests aus dem Vergleich beobachteter und unter der Nullhypothese erwarteter Ereignisse herzuleiten, den Test durchzuführen und sein Ergebnis zu interpretieren. 5. das Cox-Proportional-Hazards-Modell zu erklären, in R zu implementieren und Hazard Ratios fachgerecht zu interpretieren. 6. die zentralen Modellannahmen des Cox-Modells mithilfe von Schoenfeld- und Martingale-Residuen zu prüfen und diagnostische Auffälligkeiten zu beurteilen. 7. die Modellgüte anhand der Konkordanz zu beurteilen und aus den Ergebnissen begründete Handlungsempfehlungen zur Kundenabwanderung abzuleiten. ::: ## Theoretische Grundlagen der Survival Analysis Die Survival Analysis, auch als Überlebenszeitanalyse bekannt, beschäftigt sich mit der Analyse der Zeit bis zum Eintreten eines bestimmten Ereignisses. Im Kontext der Kundenabwanderung betrachten wir die "Überlebenszeit" eines Kunden, also die Zeitspanne von der Akquisition bis zur Kündigung eines Vertrags. ### Die Überlebensfunktion (Survival Function) Die Überlebensfunktion gibt die Wahrscheinlichkeit an, dass ein Kunde länger als den Zeitpunkt $t$ "überlebt": $$S(t) = P(T > t)$$ wobei $T$ die Überlebenszeit (Zufallsvariable) ist. Die Survival Function gibt also die Wahrscheinlichkeit an, dass das interessierende Ereignis (z. B. Kündigung, Tod, ...) **nicht** bis zum Zeitpunkt $t$ eintritt. In der Praxis ist jedoch die Abschätzung dieser Wahrscheinlichkeit besonders durch sogenannte zensierte Daten schwierig, d. h. wenn beispielsweise ein Arbeitnehmer das Unternehmen aufgrund von Krankheit, Todesfall oder Wegzug ins Ausland verlässt – also aus einem Grund, der in der Studie nicht betrachtet wird. Bei medizinischen Studien wäre dies etwa der Ausstieg einer Patientin aus der Studie. Wie können wir trotz dieser "Zensierung" mancher Teilnehmender zu zuverlässigen Schätzungen gelangen? ### Zensieren In vielen praktischen Situationen (z. B. klinische Studien, Laufzeitmessungen, Kunden"lebenszeiten") liegen zensierte Daten vor. Das heißt, für manche Probanden, Kunden oder auch Bauteile hat man nur den Zeitpunkt, bis zu dem sie beobachtet wurden (und innerhalb dieser Beobachtungszeit ist das Ereignis nicht eingetreten). Man kennt also ihre exakte „Lebensdauer“ nicht – nur, dass sie wenigstens bis zum Ende der Beobachtung durchgehalten haben. Wir unterscheiden: - **Rechtszensierung (right-censoring)**: Wenn die Studie endet oder die Person aus anderen Gründen nicht mehr beobachtet wird, ohne dass das Ereignis aufgetreten ist. - **Linkszensierung (left-censoring)**: Seltenere Variante, bei der man nur weiß, dass das Ereignis bereits vor der ersten Beobachtung eingetreten sein muss. (In der Praxis weniger typisch für Kaplan-Meier, aber theoretisch möglich.) - **Intervallzensierung (interval-censoring)**: Bei Intervallzensierung weiß man nur, dass das Ereignis in einem gewissen Zeitintervall passiert sein muss, aber nicht genau, wann. Ein praktisches Beispiel wäre, dass ein Sensor die Funktion eines Bauteils zu bestimmten Zeitpunkten misst. Wenn ein Bauteil ausfällt, kann man in der Regel nur sagen, dass es zwischen zwei Messungen ausgefallen ist. Der Kaplan-Meier-Schätzer schafft es nun, diese teilweise fehlende Information geeignet zu verarbeiten, indem er alle Personen (bzw. Objekte) bis zum Zeitpunkt, an dem wir sie verlieren, berücksichtigt und danach herausnimmt. Dies sehen wir uns jetzt an. ### Kaplan-Meier-Schätzer Statt eine globale Schätzung für die gesamte Kurve auf einmal durchzuführen, baut der Kaplan-Meier-Schätzer die Überlebenswahrscheinlichkeit **schrittweise** auf. Man stellt sich vor, dass zu Beginn alle Teilnehmenden (oder Systeme) noch „am Leben“ sind. Dann verfolgt man, was zu den Zeitpunkten der Ereignisse passiert. ::: {.callout-note title="Kernidee: Kette bedingter Wahrscheinlichkeiten"} Der Kaplan-Meier-Schätzer entsteht aus einer **Kette bedingter Überlebenswahrscheinlichkeiten**.\ Für Ereigniszeitpunkte $t_1 < t_2 < \cdots < t_k$ gilt konzeptionell: $$P(T > t_k) = P(T > t_1)\, P(T > t_2 \mid T > t_1)\, \cdots \, P(T > t_k \mid T > t_{k-1}).$$ Das heißt: Wir „überleben“ bis $t_k$, indem wir (i) bis $t_1$ überleben und (ii) danach **unter der Bedingung**, dass wir bis zum vorherigen Zeitpunkt überlebt haben, auch den nächsten Zeitpunkt überleben, usw. ::: 1. **Sortieren**\ Man sortiert alle Ereigniszeiten in aufsteigender Reihenfolge: $$t_1 < t_2 < \cdots < t_k.$$ 2. **Zählen**\ Für jeden Ereigniszeitpunkt $t_i$ bestimmt man: - $d_i$: Anzahl der Ereignisse, die genau zum Zeitpunkt $t_i$ auftreten.\ - $n_i$: Anzahl der Personen (bzw. Geräte), die kurz vor $t_i$ noch beobachtet werden (also noch nicht ausgefallen und nicht zensiert sind). Dieses Set nennt man das **Risiko-Set**. 3. **Multiplikative Aktualisierung**\ Der Kaplan-Meier-Schätzer ist die Produktform: $$\hat{S}(t) \;=\; \prod_{i : t_i \leq t} \Bigl(1 - \frac{d_i}{n_i}\Bigr).$$ - **Idee:** $\bigl(1 - d_i/n_i\bigr)$ ist die geschätzte **bedingte** Wahrscheinlichkeit, den Zeitpunkt $t_i$ zu „überleben“, **gegeben**, dass man kurz davor noch „at risk“ war.\ - Durch die Multiplikation all dieser Faktoren für alle Zeitpunkte $t_i$ bis $t$ erhält man die gesamte Überlebenswahrscheinlichkeit bis zu $t$. 4. **Schrittfunktion**\ Graphisch ergibt sich eine **Treppenfunktion**, die nur dann einen Sprung nach unten macht, wenn tatsächlich Ereignisse auftreten. Zwischen den Ereignissen bleibt sie konstant. Wir hatten die zensierten Daten als besondere Herausforderung definiert. Wie genau wird die Zensur nun beim Kaplan-Meier-Schätzer berücksichtigt? Im Schätzer werden zensierte Datenwerte durch $n_i$ berücksichtigt, d. h. in $\hat{S}(t)$ finden nur solche Kunden, Probanden, Bauteile etc. Berücksichtigung, die zum Zeitpunkt $t_i$ "at risk" waren. Damit fallen alle heraus, die z. B. aus der Studie aussteigen. $$\hat{S}(t) \;=\; \prod_{i : t_i \leq t} \Bigl(1 - \frac{d_i}{n_i}\Bigr).$$ Graphisch sehen wir dann Überlebenskurven, wie beispielsweise diese: ![](images/tier_comparison-1.png) Alle Überlebenskurven fallen ab, aber wir können mit dem Kaplan-Meier-Schätzer Unterschiede zwischen Gruppen sehen. Genau für diese Unterschiede interessieren wir uns bei der Überlebenszeitanalyse. Ob jedoch die Unterschiede statistisch signifikant sind, müssen wir im nächsten Schritt untersuchen. ### Der Log-Rank-Test Die Kaplan-Meier-Kurven beschreiben das Überlebensverhalten der Gruppen. Sichtbare Abstände zwischen den Kurven können jedoch auch durch Zufall entstehen. Der Log-Rank-Test beantwortet deshalb die Frage: > Sind die über den gesamten Beobachtungszeitraum auftretenden Unterschiede grösser, als wir sie bei gleichem Ereignisrisiko in allen Gruppen erwarten würden? #### Nullhypothese ::: {.callout-note title="Nullhypothese des Log-Rank-Tests"} Unter der Nullhypothese besitzen alle Gruppen zu jedem Zeitpunkt dieselbe Hazardfunktion: $$H_0:\quad h_1(t)=h_2(t)=\dots=h_g(t) \qquad \text{für alle }t\geq 0.$$ Bei identischem Ausgangswert $S_j(0)=1$ ist dies gleichbedeutend mit identischen Überlebensfunktionen: $$S_1(t)=S_2(t)=\dots=S_g(t).$$ ::: Der Test betrachtet ausschliesslich Zeitpunkte, an denen mindestens ein Ereignis beobachtet wird. Jeder Ereigniszeitpunkt ist dabei ein kleines bedingtes Zufallsexperiment: **Aus welcher Gruppe stammen die Ereignisse, gegeben das aktuelle Risiko-Set?** #### Das Zufallsmodell an einem Ereigniszeitpunkt Wir betrachten zunächst zwei Gruppen und einen Ereigniszeitpunkt $d_k$. Unmittelbar vor $d_k$ seien - $R_k$ das Risiko-Set, - $r_k=|R_k|$ die Gesamtzahl der Personen im Risiko-Set, - $r_{1k}$ und $r_{2k}$ die entsprechenden Anzahlen in Gruppe 1 und 2, - $q_k$ die Gesamtzahl der Ereignisse bei $d_k$, - $O_{1k}$ die tatsächlich beobachtete Anzahl dieser Ereignisse in Gruppe 1. $O_{1k}$ ist die Zufallsvariable, deren Erwartungswert wir benötigen. Unter $H_0$ besitzen alle Personen im Risiko-Set dieselbe momentane Ereignisrate. Bedingt auf $R_k$ und darauf, dass bei $d_k$ insgesamt $q_k$ Ereignisse auftreten, sind die $q_k$ betroffenen Personen daher wie eine Auswahl **ohne Zurücklegen** aus den $r_k$ Personen des Risiko-Sets. Folglich gilt $$O_{1k}\mid R_k,q_k \sim \operatorname{Hypergeom}(r_k,r_{1k},q_k).$$ #### Herleitung des Erwartungswerts Für jede Person $i$ aus Gruppe 1, die sich in $R_k$ befindet, definieren wir den Indikator $$I_{ik}= \begin{cases} 1, & \text{wenn Person }i\text{ zu den }q_k\text{ Ereignissen gehört},\\ 0, & \text{sonst}. \end{cases}$$ Dann ist $$O_{1k}=\sum_{i\in R_k\cap G_1} I_{ik}.$$ Unter $H_0$ hat jede der $r_k$ Personen dieselbe bedingte Wahrscheinlichkeit, zu den $q_k$ Ereignissen zu gehören: $$P(I_{ik}=1\mid R_k,q_k)=\frac{q_k}{r_k}.$$ Mit der Linearität des Erwartungswerts folgt $$\begin{aligned} E(O_{1k}\mid R_k,q_k) &=E\!\left(\sum_{i\in R_k\cap G_1}I_{ik}\,\middle|\,R_k,q_k\right)\\ &=\sum_{i\in R_k\cap G_1}E(I_{ik}\mid R_k,q_k)\\ &=\sum_{i\in R_k\cap G_1}\frac{q_k}{r_k}\\ &=r_{1k}\frac{q_k}{r_k}. \end{aligned}$$ Damit erhalten wir $$\boxed{E_{1k}=\frac{r_{1k}}{r_k}\,q_k.}$$ Der oft verwendete Satz „Die Ereignisse verteilen sich proportional zur Zusammensetzung des Risiko-Sets“ ist somit keine zusätzliche Regel, sondern die Interpretation dieses bedingten Erwartungswerts. #### Varianz am Ereigniszeitpunkt Da $O_{1k}$ hypergeometrisch verteilt ist, gilt $$\begin{aligned} \operatorname{Var}(O_{1k}\mid R_k,q_k) &=q_k\frac{r_{1k}}{r_k} \left(1-\frac{r_{1k}}{r_k}\right) \frac{r_k-q_k}{r_k-1}\\[2mm] &=\frac{r_{1k}r_{2k}q_k(r_k-q_k)} {r_k^2(r_k-1)}. \end{aligned}$$ Der letzte Faktor ist die Endlichkeitskorrektur der Ziehung ohne Zurücklegen. Bei nur einem Ereignis ($q_k=1$) reduziert sich die Situation auf die Frage, ob dieses eine Ereignis aus Gruppe 1 oder aus Gruppe 2 stammt. #### Aggregation über alle Ereigniszeitpunkte An jedem Ereigniszeitpunkt erhalten wir eine Abweichung $$O_{1k}-E_{1k}.$$ Positive und negative Abweichungen werden zunächst **mit ihrem Vorzeichen** über die Zeit aufsummiert: $$U=\sum_{k=1}^{K}(O_{1k}-E_{1k}).$$ Unter $H_0$ erwarten wir $U\approx 0$. Die zugehörige Varianz ist $$V=\sum_{k=1}^{K}\operatorname{Var}(O_{1k}\mid R_k,q_k).$$ Die standardisierte Log-Rank-Teststatistik für zwei Gruppen lautet deshalb $$Z=\frac{U}{\sqrt V},$$ beziehungsweise in der üblichen Chi-Quadrat-Form $$\boxed{\chi^2=\frac{U^2}{V}\sim\chi^2_1}$$ asymptotisch unter $H_0$. ::: {.callout-important title="Warum nicht jede Abweichung einzeln quadrieren?"} Der Log-Rank-Test sucht nach einer **über die Zeit konsistenten Richtung** der Abweichungen. Würden wir $(O_{1k}-E_{1k})^2$ für jeden Zeitpunkt separat summieren, gingen die Vorzeichen verloren. Auch zufällige Wechsel zwischen „zu vielen“ und „zu wenigen“ Ereignissen würden sich dann künstlich verstärken. Beim Log-Rank-Test wird deshalb zuerst die vorzeichenbehaftete Gesamtabweichung $U$ gebildet und erst danach quadriert. ::: Ein kleiner p-Wert spricht gegen die Nullhypothese gleicher Hazard- beziehungsweise Überlebensfunktionen. Der Test zeigt jedoch nicht, **wann** sich die Gruppen unterscheiden oder wie gross der Unterschied zu einem bestimmten Zeitpunkt ist. Dafür betrachten wir zusätzlich die Kaplan-Meier-Kurven und geeignete Effektmasse. ::: {.callout-tip title="Für Interessierte: mehr als zwei Gruppen" collapse="true"} Für $g$ Gruppen bilden wir einen Vektor der aufsummierten Abweichungen $$\mathbf U=(U_1,\dots,U_{g-1})^\top$$ und die zugehörige Kovarianzmatrix $\mathbf V$. Eine Gruppe ist redundant, weil sich die Abweichungen über alle Gruppen zu null summieren. Die Teststatistik lautet $$\chi^2=\mathbf U^\top\mathbf V^{-1}\mathbf U \sim\chi^2_{g-1}.$$ Damit wird zugleich verständlich, weshalb der Test bei $g$ Gruppen $g-1$ Freiheitsgrade besitzt. ::: Schliessen wir das Theoriekapitel mit zwei Verständnisfragen: ::: {.callout-note title="Verständnisfrage 1: Kaplan–Meier-Kurve – was wird eigentlich geschätzt?" collapse="false"} Eine Kaplan–Meier-Kurve fällt nur an Zeitpunkten, an denen **Ereignisse** auftreten. Zensierte Beobachtungen führen zwar zu einer kleineren Risikogruppe, aber nicht zu einem „Sprung“ der Kurve. **Frage:** Warum ist das so? Formuliere die Idee in 2–3 Sätzen (Stichworte: Risk-Set, bedingte Überlebenswahrscheinlichkeit, Zensierung). <details> <summary><strong>Antwort anzeigen</strong></summary> Die Kaplan–Meier-Schätzung basiert auf dem Produkt von **bedingten Überlebenswahrscheinlichkeiten**: Zu jedem Ereigniszeitpunkt $d_k$ wird geschätzt, welcher Anteil der Personen im Risk-Set das Ereignis *nicht* erleidet.\ Zensierung bedeutet: Wir wissen nur, dass die Person **bis zu diesem Zeitpunkt** überlebt hat, danach aber nicht mehr — daher wird sie ab diesem Zeitpunkt aus dem Risk-Set entfernt. Es passiert aber **kein Ereignis**, also gibt es keinen Grund für einen Sprung in der Überlebensfunktion. Sprünge entstehen nur durch beobachtete Ereignisse. </details> ::: ::: {.callout-note title="Verständnisfrage 2: Log-Rank-Test – Erwartung wovon?" collapse="false"} Am Ereigniszeitpunkt $d_k$ wird $$E_{1k}=\frac{r_{1k}}{r_k}\,q_k$$ berechnet. **Frage:** Welche Zufallsvariable besitzt diesen Erwartungswert, auf welche Informationen wird bedingt und welche Annahme der Nullhypothese begründet die zugrunde liegende Verteilung? <details> <summary><strong>Antwort anzeigen</strong></summary> Die Zufallsvariable $O_{1k}$ zählt, wie viele der insgesamt $q_k$ Ereignisse bei $d_k$ aus Gruppe 1 stammen. Wir bedingen auf das aktuelle Risiko-Set $R_k$ und auf die Gesamtzahl $q_k$ der Ereignisse. Unter $H_0$ haben alle Personen im Risiko-Set dieselbe Hazard. Deshalb entspricht die Gruppenzuordnung der $q_k$ Ereignisse einer Auswahl ohne Zurücklegen, also einer hypergeometrischen Verteilung. Deren Erwartungswert ist $$E(O_{1k}\mid R_k,q_k)=q_k\frac{r_{1k}}{r_k}.$$ </details> ::: Damit sind wir nun bereit, uns ein konkretes Beispiel anzusehen. ## Praktisches Beispiel implementieren Wir simulieren Daten für die Abonnementlaufzeit von Kunden eines Telekommunikationsanbieters. Die Ereigniszeiten werden dabei direkt aus einem proportionalen Hazard-Modell erzeugt. Dadurch sind Ereigniszeit, Kovariaten und Zensierung statistisch konsistent miteinander verknüpft. Die Kunden unterscheiden sich im Vertragstyp, im Serviceplan und in den monatlichen Kosten. ```{r} #| label: generate_data set.seed(123) n <- 1000 # Kundenmerkmale data <- data.frame( customer_id = 1:n, contract_type = sample( c("Monthly", "Annual"), n, replace = TRUE, prob = c(0.7, 0.3) ), service_plan = sample( c("Basic", "Standard", "Premium"), n, replace = TRUE ), monthly_charges = round(runif(n, 30, 150), 2) ) # Linearer Prädiktor des proportionalen Hazard-Modells charges_z <- as.numeric(scale(data$monthly_charges)) linear_predictor <- with(data, log(1.7) * (contract_type == "Monthly") + log(0.95) * (service_plan == "Standard") + log(0.65) * (service_plan == "Premium") + log(1.25) * charges_z ) # Ereigniszeiten aus einem PH-Modell mit konstanter Baseline-Hazard baseline_hazard <- 0.035 event_time <- rexp( n, rate = baseline_hazard * exp(linear_predictor) ) # Unabhängige Zensierung und administratives Studienende nach 36 Monaten censor_time <- runif(n, min = 12, max = 36) administrative_end <- 36 observed_time <- pmin(event_time, censor_time, administrative_end) data$churn <- as.integer( event_time <= censor_time & event_time <= administrative_end ) data$tenure <- pmax(round(observed_time, 1), 0.1) ``` ### Kaplan-Meier-Analyse ```{r} #| label: kaplan_meier # Survival-Objekt surv_obj <- Surv(time = data$tenure, event = data$churn) # Gesamte Überlebenskurve km_fit <- survfit(surv_obj ~ 1) # Stratifiziert nach Vertragstyp km_contract <- survfit(surv_obj ~ contract_type, data = data) # Plots ggsurvplot( km_fit, data = data, title = "Gesamte Überlebenskurve", xlab = "Zeit (Monate)", ylab = "Überlebenswahrscheinlichkeit", conf.int = TRUE, risk.table = TRUE ) ggsurvplot( km_contract, data = data, title = "Überlebenskurven nach Vertragstyp", xlab = "Zeit (Monate)", ylab = "Überlebenswahrscheinlichkeit", conf.int = TRUE, pval = TRUE ) ``` Die erste Überlebenszeitkurve (ohne Strata) umfasst alle Kunden des Unternehmens. Beim zweiten Plot wurden zwei Gruppen gebildet, je nach monatlicher oder jährlicher Laufzeit. Wir beobachten bis ca. 12 Monate kaum Unterschiede, jedoch erscheint danach die Anzahl Kündigungen bei Kunden mit monatlicher Laufzeit (und damit mehr Flexibilität bei der Kündigung) höher auszufallen als bei denen mit jährlicher Frist. Ob die Unterschiede jedoch statistisch signifikant sind, muss jetzt der Log-Rank-Test beantworten. ### Log-Rank-Test ```{r} #| label: logrank # Log-Rank-Test für den Vertragstyp log_rank_test <- survdiff(surv_obj ~ contract_type, data = data) print(log_rank_test) # survdiff() gibt die Teststatistik, aber nicht direkt den p-Wert aus log_rank_df <- length(log_rank_test$n) - 1 p_logrank <- pchisq(log_rank_test$chisq, df = log_rank_df, lower.tail = FALSE) p_logrank ``` Der berechnete p-Wert beträgt `r format.pval(p_logrank, digits = 3)`. Wir vergleichen ihn mit dem vorab festgelegten Signifikanzniveau $\alpha$. Ist $p<\alpha$, verwerfen wir die Nullhypothese gleicher Hazardfunktionen. Die Kaplan-Meier-Kurven benötigen wir trotzdem weiterhin: Der Test verdichtet den gesamten Verlauf auf eine Testentscheidung, während die Kurven zeigen, in welcher Richtung und in welchen Zeitbereichen die Unterschiede auftreten. Schauen wir uns die Kaplan-Meier-Plots nochmals genauer an. Neben den Unterschieden der abfallenden Treppenfunktionen sind hier noch die Bänder der Konfidenzintervalle farbig eingezeichnet. Diese werden mit zunehmender Zeit grösser, da die Unsicherheit der Schätzung zunimmt. Dies vor allem, weil die Anzahl der Teilnehmenden immer kleiner wird, d. h. wir beziehen Informationen aus einer immer kleiner werdenden Gruppe. Dies tritt ggf. auch schon früher auf, wenn wir z. B. eine kleine Gruppe mit einer großen Gruppe vergleichen. ::: {.callout-note title="Exkurs: Bootstrap bei zensierten Daten"} Bei kleinen Risiko-Sets werden Kaplan–Meier-Schätzungen und asymptotische Konfidenzintervalle zunehmend unsicher. Bootstrapverfahren können bei rechtszensierten Daten zur Schätzung von Standardfehlern oder zur Konstruktion von Konfidenzbändern eingesetzt werden. Sie ersetzen jedoch keine ausreichende Zahl beobachteter Ereignisse und sind deshalb nicht automatisch eine Lösung für sehr kleine Stichproben. Zudem muss das Resampling die Zensierungsstruktur angemessen berücksichtigen (Efron, 1981; Akritas, 1986). Dieser Ansatz geht über die Zielsetzung des Tutorials hinaus. ::: ## Fazit: Was lernen wir? Mit dem Log-Rank-Test und den Kaplan-Meier-Schätzern können wir paarweise die Überlebenszeitkurven gut vergleichen und Unterschiede evaluieren. Jedoch ist das eher explorativer Natur. Wenn wir mehrere Variablen haben, können wir auf diese Art kein komplettes Modell zusammenstellen und evaluieren, wie wir es z. B. von der multiplen linearen Regression oder der logistischen Regression gewohnt sind. Um dies zu können, wurde das Cox-Regressionsmodell (oder auch Cox Proportional Hazards Model) entwickelt, welches wir uns im nächsten Abschnitt ansehen. ## Cox-Regressionsmodell bzw. Cox Proportional Hazards Model In der Überlebenszeitanalyse interessieren wir uns für den Einfluss unserer Variablen oder auch Kovariaten auf die Überlebenszeit eines Objekts. Beim Kündigungsverhalten unserer Abokunden fragen wir uns, wie groß der Einfluss ("Hazard") der erhobenen Variablen, wie die Wahl des Abotyps oder die Höhe der monatlichen Kosten ist. Wie bei der multiplen linearen Regression möchten wir wissen, ob wir die einflussreichsten Variablen/Kovariaten identifizieren können, um ein möglichst einfaches Modell zu entwickeln. Das Ziel der Überlebenszeitanalyse ist nicht nur eine oder mehrere Kurve(n) zu zeichnen, sondern eine **Analysefrage** zu beantworten: > Welche Merkmale erhöhen oder senken die Kündigungsrate — und um wie viel? Dafür brauchen wir eine Grösse, die drei Dinge gleichzeitig kann: 1. Sie funktioniert auch mit **Zensierung** (Kunden, die am Ende der Beobachtung noch nicht gekündigt haben). 2. Sie ist eine **Momentaufnahme**: Wir wollen das Risiko *jetzt* vergleichen (nicht nur „irgendwann“). 3. Sie erlaubt einen **fairen Vergleich** zwischen Kunden, die unterschiedlich lange schon dabei sind (tenure). Einige werden jetzt einwenden, dass die Momentaufnahme mich doch gar nicht interessiert um am Ende eine Businessentscheidung zu treffen und z.B. das Abo-Modell ändere. Das ist natürlich völlig korrekt, wenn ich aber einen grösseren Zeitraum betrachte -"wie hoch ist das Risiko, dass ein bestimmter Kunde innerhalb der nächsten 12 Monate kündigt"- vermische ich den Einfluss mehrerer Effekte (z.B. Saisonalität, Länge der Kundenbeziehung). In diesem Zeitraum habe ich alle möglichen Kunden die kündigen. Kann ich jemanden, der nach 4 Monaten kündigt mit jemanden der nach 10 Monaten kündigt, vergleichen? Wahrscheinlich nicht. Betrachte ich es jedoch **zeitlokal** vergleiche ich nicht Äpfel mit Birnen, sondern Kunden, die alle zu diesem Zeitpunkt *at risk* waren. Genau das leistet die **Hazard** oder **Hazard-Funktion**. ##### Was misst die Hazard? Anschaulich: *„Unter den Kunden, die bis heute noch nicht gekündigt haben: Wie gross ist das momentane Kündigungsrisiko in der nächsten kurzen Zeitspanne?“* Damit wird das Risiko **vergleichbar**: Wir vergleichen Kunden *unter derselben Bedingung* („noch im Bestand“). ### Die Hazard-Funktion Die Hazard-Funktion beschreibt das **momentane** Risiko einer Kündigung unter der Bedingung $T > t$: $$\begin{aligned} h(t) &= \lim_{\Delta t \to 0} \frac{P(t < T \le t + \Delta t \mid T > t)}{\Delta t} \\ &= \lim_{\Delta t \to 0} \frac{P(t < T \le t + \Delta t \ \cap\ T > t)}{\Delta t \cdot P(T>t)} \\ &= \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t)}{\Delta t }\cdot \frac{1}{P(T>t)} \quad (\text{da aus } t \le T < t+\Delta t \text{ automatisch } T>t \text{ folgt})\\ &= f(t) \cdot \frac{1}{S(t)} \quad (\text{da } S(t)=P(T>t)). \end{aligned}$$ Daraus folgt dann, dass wir $f(t)$ aus $h(t)$ und $S(t)$ berechnen können: $f(t)=h(t)\cdot S(t)$. In der Praxis kann man zwar öfters auch die vereinfachende Annahme machen, dass $h(t) = \text{const.}$ gilt, wir wollen aber gleich zum flexibleren *proportional hazard model* weitergehen. ## Das Cox Proportional Hazards Model ### Die zentrale Vergleichsgröße: Hazard Ratio In der Praxis wollen wir selten das *absolute* Kündigungsrisiko zu einem Zeitpunkt $t$ wissen. Viel wichtiger ist der **Vergleich**: ::: {.callout-tip collapse="false" appearance="minimal"} Hat Person 1 zum Zeitpunkt $t$ ein höheres oder niedrigeres Kündigungsrisiko als Person 2 – und um wie viel? ::: Dazu betrachten wir die **Hazard Ratio (HR)** zweier Personen mit Kovariaten $x^{(1)}$ und $x^{(2)}$: $$\mathrm{HR}(t) = \frac{h\!\left(t\mid x^{(1)}\right)}{h\!\left(t\mid x^{(2)}\right)}.$$ Interpretation: - $\mathrm{HR}(t)=1$: gleiches **momentanes** Risiko, - $\mathrm{HR}(t)=2$: Person 1 hat das doppelte momentane Risiko wie Person 2, - $\mathrm{HR}(t)=0.5$: Person 1 hat das halbe momentane Risiko. Der Vorteil: Die HR ist eine **relative** Grösse. Sie bleibt auch dann interpretierbar, wenn das Basisrisiko über die Zeit stark schwankt (Saisonalität ("alle Arten von Fitnessstudio-Effekt im Januar"), Servicegebühren, die nur zwei Mal im Jahr in Rechnung gestellt werden etc. ). #### Ein Modell, das das Arbeiten mit Hazard Ratios einfach macht: Der Separationsansatz Damit wir den Einfluss von Kovariaten systematisch beschreiben können, wählen wir ein Modell, das die Hazard in einen **Zeitanteil** und einen **Kovariatenanteil** trennt: $$h(t\mid x)=h_0(t)\,g(x).$$ - $h_0(t)$ ist das zeitabhängige **Basisrisiko** (wie hoch ist das Kündigungsrisiko grundsätzlich über die Zeit?), - $g(x)$ beschreibt den **relativen Einfluss** der **Kovariaten**. Der entscheidende Punkt: Für zwei Personen kürzt sich das Basisrisiko in der Hazard Ratio heraus: $$\frac{h(t\mid x^{(1)})}{h(t\mid x^{(2)})} = \frac{h_0(t)g(x^{(1)})}{h_0(t)g(x^{(2)})} = \frac{g(x^{(1)})}{g(x^{(2)})}$$ Damit hängt die Hazard Ratio **nicht mehr von** $t$ ab — genau das ist die Proportional-Hazards-Idee. #### Cox Proportional Hazards: Wahl von $g(x)$ Wir wählen eine **stets positive** und gut interpretierbare Form: $$g(x)=\exp(\boldsymbol{\beta^\top x}),$$ Anm.: $\boldsymbol{\beta^T} \boldsymbol{x}$ ist ein sog. Skalarprodukt aus dem Koeffizientenvektor $\boldsymbol{\beta}$ und dem Vektor der Kovariaten $\boldsymbol{x}$ und die kompakte Schreibweise von: $$\boldsymbol{\beta}^{\top}\boldsymbol{x} = \begin{pmatrix} \beta_1 & \beta_2 & \cdots & \beta_p \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_p \end{pmatrix} = \beta_1x_1+\beta_2x_2+\cdots+\beta_px_p.$$ und erhalten das Cox-Modell $$h(t|x)=h_0(t)\exp(\beta^\top x).$$ Dann gilt für die Hazard Ratio: $$\frac{h(t|x^{(1)})}{h(t|x^{(2)})} = \exp\!\left(\beta^\top(x^{(1)}-x^{(2)})\right),$$ und insbesondere für eine Erhöhung von $x_j$ um 1 (bei konstanten übrigen Variablen): $$\mathrm{HR}=\exp(\beta_j).$$ #### Cox-Modell und Proportional Hazards – anschaulich am Kundenbeispiel In unseren Daten beschreiben die Kovariaten $\boldsymbol{x}$ Eigenschaften eines Kunden: - `contract_type` (Monthly vs. Annual) - `service_plan` (Basic / Standard / Premium) - `monthly_charges` (Monatskosten) ::: {.callout-tip collapse="false" appearance="minimal"} Die zentrale Frage lautet: **Welche dieser Kovariaten erhöhen oder senken die Kündigungsrate?** ::: ##### Was bedeutet Proportional Hazards konkret? Stellen wir uns zwei Kunden vor, die sich nur in **einer** Kovariaten unterscheiden: - Kunde A: `contract_type = Monthly` - Kunde B: `contract_type = Annual` Wenn „Monthly“ grundsätzlich riskanter ist, dann erwarten wir: $$\mathrm{HR}(t)=\frac{h(t|\text{Monthly})}{h(t|\text{Annual})} > 1.$$ **Proportional Hazards** bedeutet jetzt etwas sehr Spezifisches: > Dieses Verhältnis (die Hazard Ratio) ist **über die Zeit konstant**. Also: Der `Monthly` Kunde besitzt z. B. zu jedem Zeitpunkt ein 1.8-mal so hohes Risiko wie der `Annual`-Kunde, unabhängig davon, ob wir am Anfang der Kundenbeziehung sind oder nach vielen Monaten. ::: {.callout-important collapse="false"} Das Basisrisiko darf sich trotzdem über die Zeit ändern (z. B. viele Kündigungen zu Beginn, später weniger). Proportional heißt nur: **die Kurven sind vertikal skaliert**, nicht identisch! ::: #### Wie kommen wir an die Koeffizienten $\beta$? Für Kunde $i$ beobachten wir - $\tilde t_i$: Ereignis- oder Zensierungszeit, - $\delta_i\in\{0,1\}$: Ereignisindikator, - $x_i$: Kovariatenvektor. Ist $\delta_i=1$, kennen wir den Ereigniszeitpunkt. Ist $\delta_i=0$, wissen wir nur, dass die tatsächliche Ereigniszeit grösser als $\tilde t_i$ ist. ##### Vollständige Likelihood Der Likelihood-Beitrag einer Beobachtung lautet $$L_i=f(\tilde t_i\mid x_i)^{\delta_i} S(\tilde t_i\mid x_i)^{1-\delta_i}.$$ Mit $f(t\mid x)=h(t\mid x)S(t\mid x)$ erhalten wir äquivalent $$L_i=h(\tilde t_i\mid x_i)^{\delta_i} S(\tilde t_i\mid x_i).$$ Im Cox-Modell enthält diese vollständige Likelihood neben $\beta$ auch die unbekannte Funktion $h_0(t)$. Cox' entscheidende Idee besteht deshalb darin, nur die Information zu verwenden, die für den **relativen Kovariateneffekt** relevant ist. ##### Bedingte Wahrscheinlichkeit innerhalb des Risiko-Sets Betrachten wir einen Ereigniszeitpunkt $d_k$, an dem genau eine Person $i(k)$ kündigt. Das Risiko-Set $R(d_k)$ enthält alle Personen, die unmittelbar davor noch beobachtet und ereignisfrei sind. Bedingt darauf, dass bei $d_k$ ein Ereignis eintritt, ist die Wahrscheinlichkeit, dass gerade Person $i(k)$ betroffen ist, $$P\!\left(i(k)\text{ kündigt}\mid R(d_k),\text{ ein Ereignis bei }d_k\right) = \frac{h(d_k\mid x_{i(k)})} {\sum_{j\in R(d_k)}h(d_k\mid x_j)}.$$ Setzen wir $$h(t\mid x)=h_0(t)\exp(\beta^\top x)$$ ein, kürzt sich das für alle Personen gleiche Basisrisiko heraus: $$\frac{h_0(d_k)\exp(\beta^\top x_{i(k)})} {\sum_{j\in R(d_k)}h_0(d_k)\exp(\beta^\top x_j)} = \frac{\exp(\beta^\top x_{i(k)})} {\sum_{j\in R(d_k)}\exp(\beta^\top x_j)}.$$ Das Produkt über alle Ereigniszeitpunkte ergibt die **partielle Likelihood**: $$L_p(\beta)= \prod_{k=1}^{K} \frac{\exp(\beta^\top x_{i(k)})} {\sum_{j\in R(d_k)}\exp(\beta^\top x_j)}.$$ Wir maximieren ihre Logarithmusfunktion numerisch: $$\ell_p(\beta) = \sum_{k=1}^{K} \left[ \beta^\top x_{i(k)} - \log\!\left( \sum_{j\in R(d_k)}\exp(\beta^\top x_j) \right) \right].$$ ::: {.callout-note title="Was ist an der Likelihood partiell?"} Sie ist keine vollständige Wahrscheinlichkeitsbeschreibung der beobachteten Zeiten. Sie verwendet nur die bedingte Information darüber, **wer innerhalb des jeweiligen Risiko-Sets das Ereignis erlebt**. Genau dadurch kann $\beta$ geschätzt werden, ohne $h_0(t)$ parametrisch festzulegen. ::: Bei mehreren Ereignissen zum selben Zeitpunkt verwendet Software Näherungen beziehungsweise Tie-Verfahren wie Efron oder Breslow. ## Praktische Implementierung unseres Beispiels ### Cox Proportional Hazards Model ```{r} #| label: cox_model # Modell-Fit cox_model <- coxph( surv_obj ~ contract_type + service_plan + monthly_charges, data = data ) # Zusammenfassung print(summary(cox_model)) ``` ## Modelldiagnostik Wie bei anderen Regressionsmodellen müssen wir prüfen, ob die Struktur des Cox-Modells zu den Daten passt. Zwei Residuentypen beantworten dabei unterschiedliche Fragen: - **Schoenfeld-Residuen:** Bleibt der Effekt einer Kovariate über die Zeit konstant? - **Martingale-Residuen:** Ist die gewählte Funktionsform einer metrischen Kovariate angemessen? ### Schoenfeld-Residuen: Ist der Effekt zeitkonstant? Für jedes beobachtete Ereignis und jede Kovariate $j$ wird der Wert der tatsächlich kündigenden Person mit einem vom Modell erwarteten Wert im aktuellen Risiko-Set verglichen: $$r^{(S)}_{ij} = x_{ij}-\bar x_j(t_i;\hat\beta),$$ wobei $$\bar x_j(t_i;\hat\beta) = \frac{\sum_{\ell\in R(t_i)}x_{\ell j} \exp(\hat\beta^\top x_\ell)} {\sum_{\ell\in R(t_i)} \exp(\hat\beta^\top x_\ell)}.$$ Wenn die Proportional-Hazards-Annahme gilt, sollten die **skalierten** Schoenfeld-Residuen keinen systematischen Trend über die Zeit zeigen. Entscheidend ist deshalb nicht, ob die einzelnen Punkte stark streuen, sondern ob die geglättete Kurve zeitlich ansteigt, fällt oder sich systematisch krümmt. ::: {.callout-note title="Warum sehen die Punktwolken je nach Variablentyp anders aus?"} Bei einer binären Kovariate ist der beobachtete Wert $x_{ij}$ entweder 0 oder 1. Die Residuen bilden deshalb häufig zwei Bänder. Bei einer metrischen Kovariate kann $x_{ij}$ viele Werte annehmen; die Punkte verteilen sich entsprechend über einen grösseren vertikalen Bereich. Diese unterschiedliche Geometrie ist normal. Für die PH-Annahme zählt der **Zeittrend**, nicht die Form der Punktwolke an sich. ::: ```{r} #| label: ph_test ph_test <- cox.zph(cox_model) print(ph_test) plot(ph_test) ``` Ein kleiner p-Wert im Test von `cox.zph()` beziehungsweise ein deutlicher Trend im Plot spricht gegen einen zeitkonstanten Effekt. Ein nicht signifikanter Test beweist die PH-Annahme jedoch nicht; insbesondere bei kleinen Stichproben kann die Teststärke gering sein. ### Martingale-Residuen: Beobachtete minus erwartete Ereignisse #### Motivation: Was ist hier überhaupt ein Residuum? In einer linearen Regression ist die Idee unmittelbar: $$\text{Residuum}=\text{beobachteter Wert}-\text{vorhergesagter Wert}.$$ Bei Survivaldaten gibt es jedoch keinen vollständig beobachteten Zielwert für jede Person. Für eine zensierte Person kennen wir die Ereigniszeit gerade **nicht**. Trotzdem liefert jede Person zwei Informationen: 1. ob während ihrer Beobachtungszeit ein Ereignis auftrat, 2. wie lange sie dem Ereignisrisiko ausgesetzt war. Ein geeignetes Residuum muss daher die **beobachtete Zahl der Ereignisse** mit der **bis zum Beobachtungsende vom Modell erwarteten Zahl der Ereignisse** vergleichen. #### Schritt 1: Der beobachtete Ereignisprozess Für Person $i$ definieren wir den Zählprozess $$N_i(t) = \begin{cases} 0, & \text{solange bis }t\text{ kein Ereignis beobachtet wurde},\\ 1, & \text{nach dem beobachteten Ereignis}. \end{cases}$$ Bei einer zensierten Person bleibt $N_i(t)$ während der gesamten Beobachtung gleich 0. Zusätzlich bezeichnet $$Y_i(t)= \begin{cases} 1, & \text{wenn Person }i\text{ unmittelbar vor }t\text{ noch at risk ist},\\ 0, & \text{sonst} \end{cases}$$ ob die Person zu diesem Zeitpunkt überhaupt noch zum Risiko-Set gehört. #### Schritt 2: Was erwartet das Cox-Modell in einem kleinen Zeitintervall? Unter dem Cox-Modell besitzt Person $i$ während ihrer Risikozeit die momentane Ereignisrate $$\lambda_i(t) = Y_i(t)h_0(t)\exp(\beta^\top x_i).$$ Für ein sehr kleines Intervall $[t,t+dt)$ bedeutet dies näherungsweise $$E\!\left[dN_i(t)\mid\text{bisheriger Verlauf}\right] = \lambda_i(t)\,dt.$$ Die rechte Seite ist also der vom Modell erwartete kleine Zuwachs der Ereigniszahl. #### Schritt 3: Erwartete Ereignisse über die gesamte Beobachtungszeit Addieren beziehungsweise integrieren wir diese kleinen erwarteten Zuwächse bis zum Zeitpunkt $t$, erhalten wir $$A_i(t) = \int_0^t Y_i(u)h_0(u)\exp(\beta^\top x_i)\,du.$$ $A_i(t)$ wird **kumulative Intensität** oder **Kompensator** genannt. Inhaltlich ist es die bis $t$ aufgelaufene erwartete Ereigniszahl. Für zeitlich konstante Kovariaten und den beobachteten Endzeitpunkt $\tilde t_i$ gilt $$A_i(\tilde t_i) = H_0(\tilde t_i)\exp(\beta^\top x_i),$$ wobei $H_0(t)=\int_0^t h_0(u)\,du$ die kumulative Baseline-Hazard ist. #### Schritt 4: Beobachtet minus erwartet Nun können wir ein Residuum bilden: $$M_i(t)=N_i(t)-A_i(t).$$ Für ein kleines zukünftiges Intervall gilt bei einem korrekt spezifizierten Modell $$E\!\left[dM_i(t)\mid\text{bisheriger Verlauf}\right] = E[dN_i(t)\mid\text{bisheriger Verlauf}]-\lambda_i(t)dt =0.$$ Genau diese Eigenschaft gibt dem Prozess seinen Namen: $M_i(t)$ ist ein **Martingal**. Anschaulich besitzt die Differenz „beobachtet minus erwartet“ keinen anhand der bisherigen Information vorhersagbaren systematischen Drift. Sie kann zufällig steigen oder fallen, aber das Modell sollte keine verbleibende Richtung vorhersagen können. Am Ende der Beobachtung ist $N_i(\tilde t_i)=\delta_i$. Nach Einsetzen der geschätzten Modellgrössen erhalten wir das Martingale-Residuum $$\boxed{ \hat M_i = \delta_i- \hat H_0(\tilde t_i)\exp(\hat\beta^\top x_i) }$$ oder kurz $$\hat M_i = \text{beobachtete Ereigniszahl} - \text{erwartete Ereigniszahl}.$$ #### Interpretation eines einzelnen Martingale-Residuums - $\hat M_i>0$: Es wurde ein Ereignis beobachtet, obwohl das Modell bis zu diesem Zeitpunkt weniger als ein Ereignis erwartet hatte. Anschaulich trat das Ereignis früher oder bei geringerem vorhergesagtem kumulativem Risiko auf als vom Modell nahegelegt. - $\hat M_i\approx 0$: Beobachtete und erwartete Ereigniszahl stimmen ungefähr überein. - $\hat M_i<0$: Das Modell hat bis zum Beobachtungsende mehr Ereignisrisiko akkumuliert, als tatsächlich beobachtet wurde. Dies tritt besonders bei zensierten Personen mit hohem vorhergesagtem Risiko oder bei später als erwartet eintretenden Ereignissen auf. ::: {.callout-important title="Martingale-Residuen sind nicht symmetrisch"} Für eine Person kann höchstens ein Ereignis beobachtet werden. Deshalb gilt $\hat M_i\leq 1$. Nach unten sind die Residuen dagegen nicht entsprechend begrenzt. Zensierte Beobachtungen haben stets $\hat M_i=-\hat A_i(\tilde t_i)\leq0$. Die Residuen sind daher typischerweise schief verteilt und nicht wie gewöhnliche normalverteilte Regressionsfehler zu lesen. ::: #### Wie zeigen Martingale-Residuen eine falsche Funktionsform? Nehmen wir an, der wahre Einfluss der monatlichen Kosten sei gekrümmt, wir modellierten ihn aber als linearen Term $\beta\,\texttt{monthly\_charges}$. Dann bliebe ein Teil dieses systematischen Zusammenhangs im Residuum zurück. Eine anschauliche Diagnose besteht darin, zunächst ein Cox-Modell mit den übrigen Kovariaten, aber **ohne** die zu untersuchende metrische Variable zu schätzen. Anschliessend tragen wir die Martingale-Residuen gegen diese Variable auf: - Eine ungefähr geradlinige geglättete Beziehung spricht für einen linearen Term im Log-Hazard. - Eine deutliche Krümmung spricht für eine Transformation, einen polynomialen Term oder einen Spline. - Einzelne extreme Punkte können auf ungewöhnliche Beobachtungen hinweisen, sind wegen der asymmetrischen Residuenverteilung aber vorsichtig zu interpretieren. ```{r} #| label: martingale_residuals # Modell ohne monthly_charges: Die Residuen sollen die noch fehlende # Funktionsform dieser Kovariate sichtbar machen. cox_without_charges <- coxph( surv_obj ~ contract_type + service_plan, data = data ) martingale_plot_data <- data.frame( monthly_charges = data$monthly_charges, martingale = residuals(cox_without_charges, type = "martingale") ) ggplot( martingale_plot_data, aes(x = monthly_charges, y = martingale) ) + geom_point(alpha = 0.3) + geom_smooth(method = "loess", se = TRUE) + labs( title = "Martingale-Residuen und Funktionsform der Monatskosten", x = "Monatliche Kosten", y = "Martingale-Residuum" ) ``` Die geglättete Kurve ist hier das zentrale Diagnoseinstrument. Sie soll nicht zwingend horizontal bei null liegen, weil `monthly_charges` im Diagnosemodell bewusst fehlt. Ihre **Form** zeigt, welche funktionale Beziehung noch in das Modell aufgenommen werden sollte. Nach Wahl einer passenden Form sollte im vollständig spezifizierten Modell keine ausgeprägte systematische Reststruktur verbleiben. ::: {.callout-note title="Schoenfeld und Martingale prüfen unterschiedliche Annahmen"} - Schoenfeld-Residuen fragen: **Ändert sich der Koeffizient über die Zeit?** - Martingale-Residuen fragen: **Ist die Form des Kovariateneffekts richtig modelliert?** Eine Kovariate kann linear im Log-Hazard wirken, aber die PH-Annahme verletzen — oder umgekehrt einen zeitkonstanten, jedoch nichtlinearen Effekt besitzen. ::: ## Modellgüte im Cox-Modell (Analogien zu $R^2$) In der Survival-Analyse gibt es **kein einziges** $R^2$ wie bei der linearen Regression. Stattdessen beurteilen wir die Güte eines Cox-Modells über **Diskrimination**, **Signifikanz** und optional **Kalibrierung/Fehlermaße**. ### Diskrimination: Konkordanz (Harrell’s C) **Was misst C?** Anteil der (vergleichbaren) Paare, in denen derjenige mit **kürzerer Überlebenszeit** auch **höheren vorhergesagten Hazard** (grösseren linearen Prädiktor) hat.\ **Skala:** 0.5 = Zufall, 1 = perfekt. Eine Konkordanz von 0.5 entspricht einer zufälligen Rangordnung, 1 einer perfekten Rangordnung. Die praktische Bedeutung von Zwischenwerten ist stark vom Anwendungskontext und vom Grad der Zensierung abhängig. ```{r} #| label: concordance # Harrell's C aus der Modellzusammenfassung summary(cox_model)$concordance # Aktuelle survival-Schnittstelle concordance(cox_model) ``` ## Interpretation und praktische Empfehlungen Zunächst bereiten wir die Koeffizienten als Hazard Ratios mit Konfidenzintervallen auf: ```{r} #| label: cox_results cox_summary <- summary(cox_model) cox_results <- data.frame( Variable = rownames(cox_summary$coefficients), Hazard_Ratio = cox_summary$conf.int[, "exp(coef)"], CI_unten = cox_summary$conf.int[, "lower .95"], CI_oben = cox_summary$conf.int[, "upper .95"], p_Wert = cox_summary$coefficients[, "Pr(>|z|)"], row.names = NULL ) knitr::kable( cox_results, digits = 3, caption = "Hazard Ratios des Cox-Modells" ) ``` Für die Interpretation gelten folgende Regeln: - Eine Hazard Ratio grösser als 1 steht für eine höhere momentane Kündigungsrate gegenüber der Referenzkategorie beziehungsweise bei einer Erhöhung der Kovariate um eine Einheit. - Eine Hazard Ratio kleiner als 1 steht für eine niedrigere momentane Kündigungsrate. - Enthält das 95%-Konfidenzintervall den Wert 1, ist der entsprechende Effekt auf dem 5%-Niveau nicht statistisch signifikant. - Bei `monthly_charges` bezieht sich die Hazard Ratio auf **eine Geldeinheit**. Für eine fachlich besser verständliche Änderung um beispielsweise 10 Einheiten kann die Hazard Ratio als $\exp(10\hat\beta)$ berechnet werden. Da die Daten aus einem proportionalen Hazard-Modell simuliert wurden, sollten die geschätzten Richtungen ungefähr der Datengenerierung entsprechen: `Monthly` besitzt gegenüber `Annual` eine höhere Hazard, `Premium` gegenüber `Basic` eine niedrigere Hazard und steigende Monatskosten erhöhen das Risiko. Abweichungen der geschätzten Zahlen von den vorgegebenen Simulationseffekten entstehen durch Stichprobenzufall und Zensierung. ### Was dürfen wir daraus praktisch ableiten? Das Modell identifiziert **statistische Zusammenhänge**, aber noch keine kausalen Wirkungen. Ein höheres Kündigungsrisiko bei monatlichen Verträgen bedeutet beispielsweise nicht automatisch, dass eine Umstellung des Vertragstyps die Kündigungen im gleichen Ausmass senken würde. Für Entscheidungen können die Ergebnisse dennoch genutzt werden, um - Kundengruppen mit erhöhtem prognostiziertem Risiko genauer zu untersuchen, - gezielte Retention-Massnahmen zu testen, - mögliche nichtlineare Preiseffekte oder zeitabhängige Vertragseffekte weiter zu modellieren, - die Prognose auf unabhängigen Daten hinsichtlich Diskrimination und Kalibrierung zu validieren. ::: {.callout-important title="Trennung von Inferenz und Prognose"} Signifikante Koeffizienten zeigen Evidenz für Zusammenhänge im Modell. Eine brauchbare individuelle Vorhersage erfordert zusätzlich eine ausreichende Diskrimination, eine gute Kalibrierung und eine externe oder zumindest resampling-basierte Validierung. ::: ## Referenzen 1. Kaplan, E. L., & Meier, P. (1958). [Nonparametric estimation from incomplete observations](https://doi.org/10.1080/01621459.1958.10501452). *Journal of the American Statistical Association*, 53(282), 457–481. 2. Mantel, N. (1966). [Evaluation of survival data and two new rank order statistics arising in its consideration](https://pubmed.ncbi.nlm.nih.gov/5910392/). *Cancer Chemotherapy Reports*, 50(3), 163–170. 3. Efron, B. (1981). [Censored data and the bootstrap](https://doi.org/10.1080/01621459.1981.10477650). *Journal of the American Statistical Association*, 76(374), 312–319. 4. Akritas, M. G. (1986). [Bootstrapping the Kaplan–Meier estimator](https://doi.org/10.1080/01621459.1986.10478369). *Journal of the American Statistical Association*, 81(396), 1032–1038. 5. Cox, D. R. (1972). [Regression models and life-tables](https://doi.org/10.1111/j.2517-6161.1972.tb00899.x). *Journal of the Royal Statistical Society. Series B (Methodological)*, 34(2), 187–202. 6. Andersen, P. K., & Gill, R. D. (1982). [Cox's regression model for counting processes: A large sample study](https://doi.org/10.1214/aos/1176345976). *The Annals of Statistics*, 10(4), 1100–1120. 7. Schoenfeld, D. (1982). [Partial residuals for the proportional hazards regression model](https://doi.org/10.1093/biomet/69.1.239). *Biometrika*, 69(1), 239–241. 8. Therneau, T. M., Grambsch, P. M., & Fleming, T. R. (1990). [Martingale-based residuals for survival models](https://doi.org/10.1093/biomet/77.1.147). *Biometrika*, 77(1), 147–160. 9. Grambsch, P. M., & Therneau, T. M. (1994). [Proportional hazards tests and diagnostics based on weighted residuals](https://doi.org/10.1093/biomet/81.3.515). *Biometrika*, 81(3), 515–526. 10. Klein, J. P., & Moeschberger, M. L. (2003). [*Survival Analysis: Techniques for Censored and Truncated Data*](https://doi.org/10.1007/b97377) (2nd ed.). Springer. 11. Harrell, F. E., Jr. (2015). [*Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis*](https://doi.org/10.1007/978-3-319-19425-7) (2nd ed.). Springer.

Theoretische Grundlagen der Survival Analysis

Die Überlebensfunktion (Survival Function)

Zensieren

Kaplan-Meier-Schätzer

Der Log-Rank-Test

Nullhypothese

Das Zufallsmodell an einem Ereigniszeitpunkt

Herleitung des Erwartungswerts

Varianz am Ereigniszeitpunkt

Aggregation über alle Ereigniszeitpunkte

Praktisches Beispiel implementieren

Kaplan-Meier-Analyse

Log-Rank-Test

Fazit: Was lernen wir?

Cox-Regressionsmodell bzw. Cox Proportional Hazards Model

Was misst die Hazard?

Die Hazard-Funktion

Das Cox Proportional Hazards Model

Die zentrale Vergleichsgröße: Hazard Ratio

Ein Modell, das das Arbeiten mit Hazard Ratios einfach macht: Der Separationsansatz

Cox Proportional Hazards: Wahl von \(g(x)\)

Cox-Modell und Proportional Hazards – anschaulich am Kundenbeispiel

Was bedeutet Proportional Hazards konkret?

Wie kommen wir an die Koeffizienten \(\beta\)?

Vollständige Likelihood

Bedingte Wahrscheinlichkeit innerhalb des Risiko-Sets

Praktische Implementierung unseres Beispiels

Cox Proportional Hazards Model

Modelldiagnostik

Schoenfeld-Residuen: Ist der Effekt zeitkonstant?

Martingale-Residuen: Beobachtete minus erwartete Ereignisse

Motivation: Was ist hier überhaupt ein Residuum?

Schritt 1: Der beobachtete Ereignisprozess

Schritt 2: Was erwartet das Cox-Modell in einem kleinen Zeitintervall?

Schritt 3: Erwartete Ereignisse über die gesamte Beobachtungszeit

Schritt 4: Beobachtet minus erwartet

Interpretation eines einzelnen Martingale-Residuums

Wie zeigen Martingale-Residuen eine falsche Funktionsform?

Modellgüte im Cox-Modell (Analogien zu \(R^2\))

Diskrimination: Konkordanz (Harrell’s C)

Interpretation und praktische Empfehlungen

Was dürfen wir daraus praktisch ableiten?

Referenzen