Verbesserte Sekundäranalyse von Überlebensdaten: rekonstruktion der Daten aus veröffentlichten Kaplan-Meier-Überlebenskurven

Die Kaplan-Meier (KM) -Schätzmethode

Die Kaplan-Meier (KM) -Methode wird verwendet, um die Wahrscheinlichkeit des Auftretens des Ereignisses bis zum Zeitpunkt t, SKM (t), aus individuellen Patientendaten abzuschätzen, die aus einer RCT erhalten wurden, die einer Rechtszensur unterliegt (bei der einige Patienten für die Nachsorge verloren gehen oder am Ende des Studienzeitraums ereignisfrei sind). Das Verfahren arbeitet, indem es die IPD in Form einer Reihe von r Zeitintervallen SKM(t m) zum Ereigniszeitpunkt t m zusammenfasst :

Der Kaplan-Meier-Datenrekonstruktionsalgorithmus

Erforderliche Dateneingaben

Die erste für den Algorithmus erforderliche Eingabedatendatei enthält die extrahierten x-Achsenkoordinaten, T k, und y-Achsenkoordinaten, S k, für k = 1,…, N Punkte auf der KM-Kurve. Dazu gibt es mehrere Softwarepakete, und wir haben festgestellt, dass die Software DigitizeIt (http://www.digitizeit.de/) eine gute Leistung erbracht hat. Die KM-Kurven, extrahiert aus a.pdf-Artikel, werden in die Software eingelesen, die Achsen definiert und dann wählt der Analyst per Mausklick Punkte aus, die aus der Kurve abgelesen werden sollen. Die resultierenden T k- und S k-Koordinaten werden dann in eine Textdatei exportiert. Diese Vorarbeiten müssen sorgfältig durchgeführt werden. Die Daten sollten ausreichend sein: Jeder Schritt, der in den Abbildungen zu sehen ist, sollte während der Datenextraktion erfasst worden sein. Der Standort und die Anzahl der Klicks sind daher wichtig. Die Daten sollten auch konsistent sein: die Wahrscheinlichkeit, dass das Ereignis auftritt, nimmt mit der Zeit ab, und es sollte überprüft werden, ob dies für die extrahierten Datenpunkte immer der Fall ist. Anomalien können aufgrund der schlechten Qualität der Kurve und menschlicher Fehler bei der Steuerung der Klicks auftreten. Anomalien sollten korrigiert werden, bevor der folgende Algorithmus ausgeführt wird. Die Zeiten, zu denen die gefährdeten Zahlen in der Veröffentlichung gemeldet werden, müssen in diesen Anfangsdaten enthalten sein. Als Konvention ist der erste Datenpunkt T1 = 0 und die Wahrscheinlichkeit, das Ereignis zum Zeitpunkt 0 zu erleben, ist daher S1 = 1. Jede KM-Kurve wird separat extrahiert.

Die zweite Eingabedatendatei, die für den Algorithmus erforderlich ist, enthält Informationen zu den gemeldeten Gefährdungszahlen. Die Kurve wird aufgeteilt in i = 1,.., nint-Intervallen, für jeden haben wir die gemeldete Gefährdungszahl zu Beginn dieses Intervalls, nrisk i , der Zeitpunkt, zu dem die Gefährdungszahl bereitgestellt wird, trisk i , die erste Zeilennummer der extrahierten Koordinaten für dieses Zeitintervall lower i , und die letzte Zeilennummer der extrahierten Koordinaten für dieses Zeitintervall upper i. nrisk i und trisk i stammen aus der Originalpublikation, während lower i und upper i aus der Anzahl der Klicks in jedem Intervall stammen, um die erste Eingabedatendatei zu erstellen. Für jedes i ist das untere i gleich k, wenn T k = trisk i und das obere i gleich k, wenn Tk+1= triski+1 .

Die letzten erforderlichen Eingabedaten sind die Gesamtzahl der Ereignisse, totevents.

Wir beginnen mit der Beschreibung des Algorithmus für den Fall, dass die Anzahl der gefährdeten Personen zu Beginn der Studie und mindestens zu einem anderen Zeitpunkt gemeldet wird und wenn die Gesamtzahl der Ereignisse gemeldet wird (Fall ‚alle Informationen‘). Wir zeigen dann, wie der Algorithmus angepasst werden kann, wenn die Anzahl der gefährdeten Ereignisse nur zu Beginn der Studie gemeldet wird (Fall ’no numbers at risk‘), wenn die Gesamtzahl der Ereignisse nicht gemeldet wird (Fall ’no total events‘) und wenn keines davon gemeldet wird (Fall ’neither‘).

Der Algorithmus für den Fall ‚Alle Informationen‘

Die Anzahl der zensierten Personen ist aus den gemeldeten Daten nicht verfügbar. Wir verwenden daher die gemeldeten Risikozahlen, nrisk i. , um die Anzahl der zensierten Personen in jedem Zeitintervall ungefähr zu bestimmen i. Wir können nicht das genaue Zensurmuster innerhalb jedes Intervalls identifizieren, und so sind wir gezwungen, eine Annahme zu treffen. Wir haben angenommen, dass die Zensur innerhalb jedes der Zeitintervalle mit einer konstanten Rate erfolgt, was sinnvoll erscheint, wenn das Zensurmuster nicht informativ ist (jedes Subjekt hat eine Zensurzeit, die statistisch unabhängig von ihrer Ausfallzeit ist).

Der Algorithmus besteht aus den folgenden Schritten (auch in Abbildung 3 dargestellt).

Flussdiagramm des Algorithmus (Fall ‚alle Informationen‘).

SCHRITT 1. Wenn in Intervall i keine Personen zensiert würden, wäre die zu Beginn des folgenden Intervalls gefährdete Zahl, nris k i + 1 n o c e n s o r , die zu Beginn des Intervalls i gefährdete Zahl, multipliziert mit der Wahrscheinlichkeit, dass das Ereignis in Intervall i auftritt, vorausgesetzt, dass es zu Beginn des Intervalls i am Leben ist:

auf die nächste ganze Zahl gerundet.

Unsere anfängliche Vermutung für die in Intervall i zensierte Zahl ist die Differenz zwischen der gemeldeten gefährdeten Zahl zu Beginn des Intervalls i + 1, nriski+1, und der gefährdeten Zahl ohne Zensur:

SCHRITT 2. Wir verteilen das c=1,…,nce n ^ so r i Zensierzeiten, ce n ^ t c , gleichmäßig über Intervall i:

Die Anzahl der zensierten Beobachtungen zwischen den extrahierten KM-Koordinaten k und k + 1 wird ermittelt, indem die Anzahl der geschätzten Zensorzeiten ce n ^ t c gezählt wird, die zwischen der Zeit T k und Tk + 1 liegen:

wobei I {c e n ^ t c ∈ } ein Indikator ist, der 1 zurückgibt, wenn ce n ^ t c auf dem Intervall liegt, und andernfalls 0.

SCHRITT 3. Die Anzahl der Ereignisse, d ^ k , an jeder extrahierten KM-Koordinate, k, und damit die Anzahl der Risikopatienten an der nächsten Koordinate, n ^ k + 1 , kann dann berechnet werden. Neu arrangieren Eq. 2 erhalten wir, dass d ^ k gleich der Anzahl der Risikopatienten bei der extrahierten KM-Koordinate k ist, multipliziert mit eins minus der Wahrscheinlichkeit, das Ereignis bei der extrahierten KM-Koordinate k zu erleben, geteilt durch Ŝ l a s t ( k) K M die geschätzte KM-Überlebenswahrscheinlichkeit bei der vorherigen Koordinate, bei der wir schätzen, dass ein Ereignis aufgetreten ist, zuletzt (k). Die Intervalle der KM-Schätzungen sind so ausgelegt, dass zu Beginn jedes Intervalls mindestens ein Ereignis auftritt, Dies ist jedoch nicht unbedingt für unsere extrahierten Koordinaten der Fall, Daher müssen wir die Zeit des letzten Ereignisses verfolgen:

wobei k‘ so ist, dass d ^ k ‚ >0

aber d ^ j = 0für j = k’+ 1,…, k – 1

Mit gl.2, wir haben:

Daher:

auf die nächste ganze Zahl gerundet.

Die Anzahl der Risikopatienten an jeder extrahierten Koordinate k wird dann unter Verwendung von Gl.1:

wobei wir zu Beginn des Intervalls n ^ l o w e r i = n ^ k i setzen. Dies ergibt eine geschätzte Anzahl von Risiken zu Beginn des folgenden Intervalls nrîs k i + 1 = n ^ u p p e r i + 1 .

SCHRITT 4. Wenn nrîs k i + 1 ≠nris k i + 1 dann passen wir die geschätzte Anzahl der zensierten Beobachtungen im Intervall i, ncenŝor um:

Wir wiederholen die Schritte 2-3 iterativ, bis die geschätzte und die veröffentlichte Risikozahl übereinstimmen (dh nrîs k i + 1 = nris k i + 1 ).

SCHRITT 5. Wenn i + 1 nicht das letzte Intervall ist, wiederholen wir die Schritte 1-4 für das folgende Intervall.

SCHRITT 6. In veröffentlichten RCTs wird im Allgemeinen keine Risikozahl am Ende des letzten Intervalls, nint, veröffentlicht. Wir nehmen zunächst an, dass die im letzten Intervall zensierte Anzahl gleich der vor dem letzten Intervall geschätzten Gesamtzahl der zensierten ist, ∑ i = 1 n i n t – 1 n c e n ŝ o r i , gewichtet mit der verbleibenden Zeit relativ zur bereits verstrichenen Zeit, gerundet auf die nächste ganze Zahl. Wenn diese Zahl jedoch größer war als die Anzahl der Patienten, die zu Beginn des letzten Intervalls noch gefährdet waren, wurde stattdessen diese Risikozahl gewählt. Diese Annahme ist formal in der folgenden Gleichung geschrieben:

Und wir führen Schritt 2-3 aus.

SCHRITT 7. Wir verwenden dann die gemeldete Gesamtzahl der Ereignisse, totevents. Wir berechnen die geschätzte Gesamtzahl der Ereignisse, die zu Beginn des letzten Intervalls erhalten wurden, ∑ k = 1 u p p e r n i n t – 1 d ^ k. Wenn dies größer oder gleich totevents ist, gehen wir davon aus, dass keine Ereignisse oder Zensuren mehr auftreten:

SCHRITT 8. Wenn ∑ k = 1 u p p e r n i n t – 1 d ^ k kleiner als totevents ist, passen wir die geschätzte Anzahl der zensierten Beobachtungen im Intervall nint, nce n ^ so r n i n t um die Differenz der Gesamtzahl der Ereignisse an:

Wir wiederholen dann die Schritte 2-3,8 für das letzte Intervall, nint, bis die geschätzte Gesamtzahl der Ereignisse, ∑ k = 1 u p p e r n i n t – 1 d ^ k , gleich der gemeldeten Gesamtzahl der Ereignisse ist, totevents oder bis die geschätzte Gesamtzahl der Ereignisse kleiner ist als die gemeldete Gesamtzahl der Ereignisse, aber die Gesamtzahl der Zensuren im letzten Intervall, nce n ^ so r n i n t , wird gleich Null.

Anpassungen des Algorithmus für den Fall ‚keine Zahlen gefährdet‘

In diesem Fall gibt es nur ein Intervall nint = 1. Wir gehen zunächst davon aus, dass die Gesamtzahl der zensierten gleich Null ist, und gehen dann wie in Schritt 8 vor.

Anpassungen des Algorithmus für den Fall ‚Keine Gesamtereignisse‘

In diesem Fall gehen wir wie für den Fall ‚Alle Informationen‘ vor, mit der Ausnahme, dass keine erneute Anpassung anhand der Gesamtzahl der Ereignisse vorgenommen werden kann und wir daher bei Schritt 6 anhalten.

Anpassung des Algorithmus für den ‚weder‘ -Fall

Wenn weder die Gesamtzahl der Ereignisse noch die Anzahl der gefährdeten Ereignisse über den Beginn der Studie hinaus gemeldet werden, gingen wir davon aus, dass es keine zensierten Beobachtungen gab. Dies ist eine starke Annahme, aber so stark wie jede andere Annahme, die wir über die Zensur ohne weitere Informationen machen könnten. Aufgrund des Mangels an Informationen wird eine geringere Qualität der Ergebnisse erwartet.

Erhalten der individuellen Patientendaten (IPD) aus den rekonstruierten Kaplan-Meier-Daten

Aus unseren rekonstruierten Kaplan-Meier-Parametern d ^ k, c, n k, n ^ k für jede extrahierte KM-Koordinate k = 1,…, N, können wir die IPD ableiten, die diese Daten erzeugen würde. Dieses letzte Stück Codierung ist in der Tat ziemlich einfach. Jedes Mal, wenn ein Ereignis oder eine Zensur geschätzt wird, wird die entsprechende Zeit sowie ein Ereignisindikator (eins für Ereignis und Null für Zensur) aufgezeichnet.

Bewertung der Reproduzierbarkeit und Genauigkeit

In der Validierungsübung wurden sechs Kaplan-Meier-Kurvenpaare verwendet. Diese wurden aus einer Teilmenge von Publikationen gezogen, die Teil eines Rückblicks auf Überlebenszeitanalysemethoden waren, die in wirtschaftlichen Bewertungen verwendet wurden . Wir führten eine Rekonstruktion von zweiundzwanzig Überlebenswahrscheinlichkeiten, sieben mittleren Überlebenszeiten, sechs Hazard Ratios und vier Standardfehlern der Log Hazard Ratios durch, die in diesen vier Publikationen berichtet wurden. Jeder wurde zweimal von denselben drei Beobachtern rekonstruiert. Zwei der drei Beobachter waren nicht an der Entwicklung des Algorithmus beteiligt.

Die Reproduzierbarkeit und Genauigkeit der Methode wurde für jede der 4 verschiedenen Informationsebenen bewertet (‚alle Informationen‘, ‚keine gefährdeten Zahlen‘, ‚keine Gesamtereignisse‘ und ‚weder noch‘). Um die Unterschiede zwischen den rekonstruierten Statistiken und den ursprünglichen zu bewerten, wurde die natürliche Skala für die Überlebenswahrscheinlichkeiten verwendet, während die Log-Skala für Mediane, HRs und ihre Unsicherheiten verwendet wurde. Kaplan-Meier-Kurven und Cox-HRs basierend auf rekonstruierten Daten wurden unter Verwendung der R-Routinen survfit und coxph geschätzt.

Wir haben eine Standard-Zwei-Wege-ANOVA mit wiederholten Messungen an die Unterschiede zwischen den rekonstruierten Ergebnissen und den ursprünglichen Ergebnissen angepasst, entweder auf der natürlichen oder der Log-Skala, abhängig von der betrachteten Statistik. Die Komponenten der Varianz waren Exemplar, Beobachter, Exemplar × Beobachter-Interaktion und innerhalb der Zelle Fehler. Da der p-Wert aus dem F-Ratio-Test für die Interaktion in allen Fällen über 10% lag, haben wir den Interaktionsterm mit dem zelleninternen Fehlerterm gepoolt. Der gewählte Ansatz ähnelt dem, was in technischen Anwendungen als ‚Messwiederholbarkeit und Reproduzierbarkeit‘ bezeichnet wird .

Die Reproduzierbarkeit stellt den Fehler dar, wenn ein einzelner Beobachter eine einzelne Rekonstruktion für eine bestimmte Statistik durchführt. Dies wurde als Summe des Within-Observer- und Between-Observer-Fehlers geschätzt. Die Monte-Carlo-Simulation aus dem angepassten ANOVA-Modell wurde verwendet, um die 95% -Konfidenzintervalle um die Standardabweichungen zu erhalten. Es wurde angenommen, dass die Freiheitsgrade für die Variationen innerhalb, zwischen und Ergebnis Chi-Quadrat-Verteilungen folgen. Um eine robuste Inferenz zu gewährleisten, wurden 150 000 Stichproben von Freiheitsgraden aus jeder dieser Verteilungen gezogen, d. H. Für jede Variationsquelle. Anschließend wurden die Mittelwertquadratschätzungen unter Verwendung der Summe der durch die ANOVA erhaltenen Quadrate und der durch die Simulation erhaltenen Stichprobe für jede der 150 000 Stichproben und für jede der Variationsquellen berechnet. Die entsprechenden 150 000 innerhalb, zwischen und Ergebnis Standardabweichungen wurden anschließend geschätzt, und wir extrahierten schließlich die 2.5 und 97.5 Perzentile, um die Konfidenzintervalle Schätzungen zu erhalten.

Um die Genauigkeit zu beurteilen, untersuchten wir die mittlere Differenz zwischen den rekonstruierten Statistiken und den ursprünglichen. Die resultierende mittlere Verzerrung oder der mittlere Fehler (MEAN Error, ME) spiegelt eine systematische Über- oder Unterschätzung wider. Die 95% -Konfidenzintervalle werden direkt aus der Schätzung der Standardabweichungen der ANOVA erhalten. Wir haben auch absolute Bias oder Mean Absolute Error (MAE) aufgezeichnet. Dies ignoriert die Richtung der Fehler und misst ihre Größe, was ein Maß für die absolute Genauigkeit der rekonstruierten Ergebnisse gibt. Eine Simulationsmethode wurde erneut verwendet, um die 95% -Konfidenzintervalle zu erhalten, die davon ausgingen, dass MEs normalverteilt waren. Um eine robuste Schlussfolgerung zu gewährleisten, wurden für jede Statistik 150 000 Stichproben aus der Normalverteilung mit dem beobachteten Mittelwert und der Varianz gemäß der ANOVA gezogen. Wir berechneten dann die entsprechenden 150 000 absoluten Werte dieser Zahlen und extrahierten schließlich die 2,5 und 97,5 Perzentile, um die Konfidenzintervallschätzungen zu erhalten.

Schließlich haben wir die Variation der Differenz zwischen rekonstruierter und ursprünglicher Statistik aufgezeichnet, die auf die Wahl der Beispiele zurückzuführen war, d. H. auf die 22 Überlebenswahrscheinlichkeiten, 7 Mediane, 6 Stunden und 4 Standardfehler der Log-Stunden. Dies gibt einen weiteren Hinweis auf die Genauigkeit des Verfahrens.