Überlebensanalyse mit Python 2-Kaplan-Meier-Schätzung

Klicken Sie hier für den vorherigen Artikel Mit Python 1 gelernte Überlebenszeitanalyse - Was sind Überlebenszeitdaten? https://qiita.com/Goriwaku/items/8d00696d853da73505bd

Informationen zu den diesmal verwendeten Daten finden Sie im vorherigen Artikel oder zum Herunterladen (whas500.csv) von meinem Git (https://github.com/goriwaku/survival_analysis).

Was ist eine Überlebensfunktion?

Ich werde die Überlebensfunktion mit WHAS500 erklären, die im vorherigen Artikel besprochen wurde. Für die Zeit t als kontinuierliche Variable definieren wir die Überlebenszeit als stochastische Variable als T. Zu diesem Zeitpunkt zeigt die kumulative Verteilungsfunktion von T an, dass die zufällig extrahierten Subjekte kleiner oder gleich einer bestimmten Zeit t sind und wie folgt definiert sind.

F(t)=Pr(T \leq t)

Hier kann die Überlebensfunktion durch die Wahrscheinlichkeit ausgedrückt werden, eine Überlebenszeit T zu erhalten, die größer als eine bestimmte Zeit t ist.

S(t)=Pr(T > t)

Es wird sein. Zu diesem Zeitpunkt, zu jedem Zeitpunkt t, gilt unter Berücksichtigung des Überlebens als verbleibendes Todesereignis die folgende Gleichung nach dem Axiom der Wahrscheinlichkeit.

S(t)=1-F(t)

Schauen wir uns nun diese Überlebensfunktionsschätzung an, insbesondere die Kaplan-Meier-Schätzung, bei der es sich um eine Schätzung der Überlebenszeitfunktion bei Vorhandensein von Verkürzungen auf der rechten Seite handelt.

Überlebensfunktionsschätzung

Wir werden die Überlebensfunktion tatsächlich schätzen. Betrachten Sie als Beispiel die Schätzung der Überlebensfunktion anhand der Überlebenszeit und des Überlebensstatus von 5 Personen, wie unten gezeigt.

Zum Zeitpunkt 0 sind alle 5 am Leben, also $ \ hat {S} (t) = 1,0 $. Da Subjekt Nummer 1 am 6. Tag stirbt, ist die Bedingung für das Sterben in diesem kleinen Abschnitt unter Berücksichtigung des Minutenabschnitts $ (6- \ delta, 6] $, der kurz vor dem 6. Tag beginnt und am 6. Tag endet Die geschätzte Wahrscheinlichkeit einer Bindung beträgt $ 1/5 $ und die Überlebenswahrscheinlichkeit beträgt $ 4/5 $. Personen, die zu einem bestimmten Zeitpunkt am Leben sind, sind ** Risikosätze ** mit einem Todesrisiko. Die Anzahl der Personen wird als gefährdete Anzahl ausgedrückt. Die geschätzte Überlebenszeit kann ausgedrückt werden durch ** (Wahrscheinlichkeit, bis zu diesem Punkt zu leben) x (bedingte Überlebenswahrscheinlichkeit im obigen Minutenabschnitt) **. Das Intervall $ I_i $ wird durch die Zeit in der obigen Tabelle in aufsteigender Reihenfolge geteilt (Beispiel: $ I_0 = [0, 6) $), $ n_i $ ist die zu diesem Zeitpunkt gefährdete Zahl und $ d_i $ ist zu diesem Zeitpunkt. Die Anzahl der Menschen, die gestorben sind. Zu diesem Zeitpunkt kann die geschätzte Sterblichkeitswahrscheinlichkeit als $ d_1 / n_1 $ und die Überlebenswahrscheinlichkeit als $ (n_1 - d_1) / n_1 $ in einem kleinen Intervall nahe dem Tag ausgedrückt werden, an dem ein Ereignis auftritt. Wenn daher die Überlebenswahrscheinlichkeit am 6. und 14. Tag geschätzt wird,

\hat{S}(6)=1 \times \frac{4}{5}=\frac{4}{5}\\
\hat{S}(14)=1 \times \frac{4}{5} \times \frac{3}{4} = \frac{3}{5}

Es wird sein. Im nächsten Fall wird Subjekt Nr. 3 mit einer nicht todunabhängigen Kündigung aus dem am 21. Tag festgelegten Risiko ausscheiden. In diesem Minutenintervall ist $ n_3 = 3 $ und $ d_3 = 0 $, also ist der geschätzte Wert der Überlebensfunktion

\hat{S}(21)=1 \times \frac{4}{5} \times \frac{3}{4} \times \frac{3}{3} = \frac{3}{5}

Es überrascht nicht, dass sich die Schätzungen der Überlebensfunktionen nicht ändern, da im Intervall $ I_3 $ keine Todesfälle auftreten. Das folgende Intervall kann auch verwendet werden, um den geschätzten Wert der Überlebensfunktion zu erhalten. Der geschätzte Wert der Überlebensfunktion, der mit dieser Methode erhalten wird, wird als ** Kaplan-Meier-Schätzwert ** bezeichnet.

Zeichnen wir nun die Kaplan-Meier-Schätzung als Schrittfunktion.

test_data = [[6, 1],[44, 1],[21, 0],[14, 1],[62, 1]]
test_data = sorted(test_data, key=lambda x: x[0], reverse=False)

s = 1
n = 5
pre = 0
for t, censor in test_data:
    plt.plot((pre, t), (s, s), color='blue')
    if censor == 1:
        plt.plot((t, t), (s, s * (n - 1) / n), color='blue')
        s = s * (n - 1) / n
    n -= 1
    pre = t
plt.show()

Ich habe die Tabelle als zweidimensionales Array gespeichert und sie zur Vereinfachung in aufsteigender Reihenfolge nach t sortiert. Wenn Sie ein zweidimensionales Array sortieren möchten, geben Sie an, was nach Schlüsselwortargumenten wie key = lambda x: x [0] sortiert werden soll. Wenn Sie nach Zensor sortieren möchten, können Sie x [1] eingeben. Danach wird die Überlebenswahrscheinlichkeit mit der for-Anweisung berechnet und die grafische Darstellung ist in der folgenden Abbildung dargestellt. Somit wird die Kaplan = Meier-Kurve als Sprungfunktion ausgegeben. Sie nimmt ab, wenn die Mortalität beobachtet wird, und bleibt während dieser Zeit konstant. Bei $ t = 62 $ gibt es keine Überlebenden, daher ist das Ende $ \ hat {S} (62) = 0 $. In diesem Beispiel sind nicht mehrere Todesfälle im selben t aufgetreten (dies wird als ** Gleichstand ** bezeichnet), aber selbst wenn ein Gleichstand auftritt, wird eine Zufallszahl generiert und die Reihenfolge wird zufällig zugewiesen. Sie können eine Diskussion entwickeln. Selbst wenn die Verbindungsdaten in der KM-Schätzung einheitlich behandelt werden, ist die endgültige Schätzung dieselbe, sodass grundsätzlich keine Anpassungen erforderlich sind. In der Praxis ist dies selten der Fall, aber wenn Sie eine extrem große Anzahl von Bindungen haben, sollten Sie in Betracht ziehen, anstelle einer KM-Schätzung ein diskretes Zeitmodell zu verwenden. Es sollte auch beachtet werden, dass die geschätzte Überlebenszeit danach nicht definiert werden kann, wenn die endgültige Beobachtungszeit auf der rechten Seite abgeschnitten wird.

Verallgemeinerung von Kaplan-Meier-Schätzungen

Kaplan-Meier-Schätzungen werden in der Überlebensanalyse so häufig verwendet, dass wir hier eine allgemeine Formulierung geben. Wenn t endlich ist, verliert die Sortierung nicht ihre Allgemeingültigkeit. Sortieren Sie also t in aufsteigender Reihenfolge und die binären Variablen $ c_i $, die bei $ t_i $ einem Risiko ausgesetzt sind, unabhängig davon, ob es sich um einen Sensor aufgrund des Todes bei $ t_i $ handelt oder nicht. Unter $ n_1 $ ist die Anzahl der beobachteten Todesfälle $ d_i $, die Kaplan-Meier-Schätzung der Überlebensfunktion zum Zeitpunkt t ist:

\hat{S}(t)=\prod_{t_i \leq t}\frac{n_i - d_i}{n_i}

Wenn jedoch $ t \ leq t_1 $ ist, ist $ \ hat {S} (t) = 1 $.

Kaplan-Meier-Kurvenplot - Lebenslinienausgabe

Ich habe die Kaplan-Meier-Kurve früher selbst gezeichnet, aber Python verfügt über eine Bibliothek zur Überlebenszeitanalyse namens Lebenslinien, ohne sie selbst implementieren zu müssen. In diesem Abschnitt werden wir diese Lebenslinien verwenden, um die Kaplan-Meier-Kurve des vorherigen WHAS500 zu zeichnen. Lebenslinien sind standardmäßig nicht enthalten. Wenn Sie sie also noch nie verwendet haben, führen Sie zunächst die folgenden Schritte in einem Terminal oder einer Eingabeaufforderung aus:

pip install lifelines

Hier ist der Python-Code.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lifelines as ll
from lifelines import KaplanMeierFitter


whas = pd.read_csv('whas500.csv')
kmf = KaplanMeierFitter()
kmf.fit(whas.LENFOL, event_observed=whas.FSTAT)
kmf.plot_survival_function()
print(kmf.survival_function_)
plt.show()

Geben Sie die Zeit t in das erste Anpassungsargument von KaplanMeierFitter of lifelines ein und geben Sie eine binäre Variable ein, ob das Ereignis (Tod), das Sie beobachten möchten, im Schlüsselwortargument event_observed aufgetreten ist, und lassen Sie es passen. Anschließend können Sie die KM-Kurve einfach mit der Funktion .plot_survival_fuction () des Modells zeichnen. Verwenden Sie außerdem Survival_function_, um die geschätzte Überlebensfunktion abzurufen. Es handelt sich um einen DataFrame-Pandas-Typ, der ein Paar aus Zeitpunkt und Schätzwert zurückgibt. Es ist auch möglich, die oben erwähnte kumulative Verteilungsfunktion $ F (t) $ zu erhalten, indem die Überlebensfunktion durch die kumulative Dichte ersetzt wird. Das tatsächliche Diagramm der KM-Kurve ist unten dargestellt. Die Breite von Hellblau in der Abbildung ist das Konfidenzintervall.

Damit ist die Erläuterung der Kaplan-Meier-Schätzungen und deren Umsetzung abgeschlossen. Das nächste Mal möchte ich die Kaplan-Meier-Schätzung interpretieren und die beiden Gruppen vergleichen. Ich hoffe, Sie werden weiterlesen.

Referenzen / Referenzlinks

Einführung in die Überlebenszeitanalyse Hosmer DW, Lemeshow S, May S. LIFELINES　https://lifelines.readthedocs.io/en/latest/index.html Kyoto University OCW Auditkurs der Kyoto University Graduate School of Medicine Biostatistik für klinische Forscher "Grundlagen der Überlebenszeitanalyse" https://www.youtube.com/watch?v=NmZaY2tDKSA&feature=emb_title