Enhanced secondary analysis of survival data: rekonstrukcja danych z opublikowanych krzywych przeżycia Kaplana-Meiera

metoda szacowania Kaplana-Meiera (KM)

metoda Kaplana-Meiera (KM) służy do oszacowania prawdopodobieństwa wystąpienia zdarzenia aż do czasu t, SKM(T), na podstawie indywidualnych danych pacjentów uzyskanych z RCT, który podlega cenzurze praw (w przypadku gdy niektórzy pacjenci są traceni do obserwacji lub nie występują zdarzenia na koniec okresu badania). Metoda działa poprzez podsumowanie IPD w postaci serii R przedziałów czasowych SKM (t m) w czasie zdarzenia t m :

algorytm rekonstrukcji danych Kaplana-Meiera

wymagane dane wejściowe

pierwszy plik danych wejściowych wymagany dla algorytmu zawiera wyodrębnione współrzędne osi x, T k i współrzędne osi y, S K, dla k = 1,…, N punktów na krzywej KM. Istnieje kilka pakietów oprogramowania, aby to zrobić, i odkryliśmy, że oprogramowanie DigitizeIt (http://www.digitizeit.de/) działa dobrze. Km krzywe, wydobyte zartykuł pdf, są wczytywane do oprogramowania, osie są zdefiniowane, a następnie analityk za pomocą kliknięć myszką wybiera punkty do odczytania z krzywej. Otrzymane współrzędne T k I S K są następnie eksportowane do pliku tekstowego. Ta wstępna praca musi być wykonana ostrożnie. Dane powinny być wystarczające: każdy krok widoczny na rysunkach powinien zostać uchwycony podczas ekstrakcji danych. Lokalizacja i liczba kliknięć są zatem ważne. Dane powinny być również spójne: prawdopodobieństwo wystąpienia zdarzenia maleje z czasem i należy sprawdzić, czy zawsze ma to miejsce w przypadku wyodrębnionych punktów danych. Anomalie mogą wystąpić ze względu na jakość publikacji krzywej i błąd ludzki w kontrolowaniu kliknięć. Wszelkie anomalie należy skorygować przed uruchomieniem poniższego algorytmu. W tych wstępnych danych należy uwzględnić czasy, w których liczby zagrożone są zgłaszane w publikacji. Zgodnie z konwencją, pierwszym punktem danych jest T1 = 0, a prawdopodobieństwo wystąpienia zdarzenia do czasu 0 wynosi zatem S1 = 1. Każda krzywa KM jest wyodrębniana oddzielnie.

drugi plik danych wejściowych wymagany dla algorytmu zawiera informacje o zgłoszonych zagrożonych liczbach. Krzywa dzieli się na i = 1,.., nint intervals, dla każdego mamy zgłoszoną liczbę zagrożoną na początku tego przedziału, nrisk i, czas, w którym podana jest liczba zagrożona , trisk i , numer pierwszego wiersza wyodrębnionych współrzędnych dla tego przedziału czasu niższego i oraz numer ostatniego wiersza wyodrębnionych współrzędnych dla tego przedziału czasu górnego i . nrisk i i trisk i pochodzą z oryginalnej publikacji, podczas gdy dolne i i górne pochodzą z liczby kliknięć wykonanych w każdym przedziale, w celu utworzenia pierwszego pliku danych wejściowych. Dla każdego i, dolne i jest równe k, gdy t k = trisk i i górne i jest równe k, gdy TK+1= triski+1.

ostateczne wymagane dane wejściowe to całkowita liczba zdarzeń, totevents.

zaczynamy od opisania algorytmu dla przypadku, w którym liczba zagrożona jest zgłaszana na początku badania i co najmniej w jednym innym punkcie czasowym i kiedy zgłaszana jest całkowita liczba zdarzeń (przypadek „wszystkie informacje”). Następnie pokazujemy, w jaki sposób algorytm można dostosować, gdy liczba zagrożona jest zgłaszana tylko na początku badania (przypadek „brak liczby zagrożonej”), gdy całkowita liczba zdarzeń nie jest zgłaszana (przypadek „brak całkowitych zdarzeń”) i gdy żadne z tych zdarzeń nie jest zgłaszane (przypadek „żaden”).

algorytm dla sprawy „wszystkie informacje”

liczba ocenzurowanych osób nie jest dostępna na podstawie zgłoszonych danych. W związku z tym wykorzystujemy zgłoszone liczby zagrożone, nrisk i, w celu przybliżenia liczby ocenzurowanych osób w każdym przedziale czasowym i. Nie możemy zidentyfikować dokładnego wzorca cenzurowania w każdym przedziale, więc jesteśmy zmuszeni do założenia. Przyjęliśmy, że cenzurowanie odbywa się w stałym tempie w każdym z przedziałów czasowych, co wydaje się rozsądne, jeśli schemat cenzurowania nie ma charakteru informacyjnego (każdy podmiot ma czas cenzurowania, który jest statystycznie niezależny od czasu jego niepowodzenia).

algorytm składa się z następujących kroków (również zilustrowanych na rysunku 3).

schemat algorytmu (przypadek 'wszystkie informacje’).

punkt 1. Najpierw tworzymy wstępne przypuszczenie dla liczby ocenzurowanej na interwale i. jeśli nie było osób ocenzurowanych na interwale i, to liczba zagrożona na początku następnego interwału, nris k i + 1 n o c e N s O r, byłaby liczbą zagrożoną na początku interwału i, pomnożoną przez prawdopodobieństwo przeżycia zdarzenia w interwale i, pod warunkiem bycia żywym na początku interwału i:

w zaokrągleniu do najbliższej liczby całkowitej.

nasze wstępne przypuszczenie dotyczące liczby ocenzurowanej w przedziale i jest różnicą między zgłoszoną liczbą zagrożoną na początku przedziału i + 1, nriski+1, A liczbą zagrożoną w przedziale bez cenzury:

punkt 2. Rozkładamy C=1,…, nce n ^ więc R I C, ce n ^ t c, równomiernie na przedział i:

liczbę ocenzurowanych obserwacji między współrzędnymi km K I k + 1 można znaleźć, licząc liczbę szacowanych czasów cenzury, ce n ^ T c, które leżą między czasem T k i Tk+1:

gdzie I { c e N ^ T c ∈} jest wskaźnikiem zwracającym 1, jeśli ce n ^ t c leży w przedziale, a 0 w przeciwnym razie.

Krok 3. Następnie można obliczyć liczbę zdarzeń, D ^ k, przy każdym pobranym koordynacie KM, k, a tym samym liczbę pacjentów zagrożonych przy następnym koordynacie, n ^ k + 1. Re-arranging Eq. 2, otrzymujemy, że d ^ k jest równe liczbie pacjentów zagrożonych przy pobranym km współrzędnych, K, pomnożonej przez jeden minus prawdopodobieństwo wystąpienia zdarzenia przy pobranym km współrzędnych, K, podzielonej przez Ŝ l A s t (k) K M szacowane prawdopodobieństwo przeżycia KM przy poprzednim współrzędnym, gdzie szacujemy, że zdarzenie miało miejsce, last(K). Przedziały oszacowań KM są zaprojektowane tak, aby co najmniej jedno zdarzenie miało miejsce na początku każdego przedziału, ale niekoniecznie jest tak w przypadku naszych wyodrębnionych współrzędnych, dlatego musimy śledzić czas ostatniego zdarzenia:

gdzie k’ jest takie, że d ^ k ’>0

ale d ^ j = 0for j = k ’ + 1,…, k-1

używając korektora.2, mamy:

:

zaokrąglone do najbliższej liczby całkowitej.

liczbę pacjentów zagrożonych każdym wyekstrahowanym koordynatem, k, uzyskuje się za pomocą Eq.1:

gdzie na początku przedziału ustawiamy n ^ L o w e r i =nris k i . Daje to szacunkową liczbę zagrożonych na początku następującego przedziału nrîs k i + 1 = n ^ u p p e R i + 1 .

KROK 4. Jeśli nrîs k i + 1 ≠nris k i + 1, to ponownie dostosowujemy szacunkową liczbę ocenzurowanych obserwacji w przedziale i, ncenŝor, przez:

powtarzamy iteracyjnie kroki 2-3, aż do oszacowania i opublikowania dopasowania liczby zagrożonej (np. nrîs k i + 1 =nris k i + 1).

krok 5. Jeśli i + 1 nie jest ostatnim interwałem, powtarzamy kroki 1-4 dla następnego interwału.

Krok 6. W opublikowanych RCTs nie ma na ogół numeru zagrożonego opublikowanego na końcu ostatniego interwału, nint. Najpierw Zakładamy, że liczba ocenzurowana w ostatnim przedziale jest równa całkowitej liczbie ocenzurowanej oszacowanej przed ostatnim przedziałem, ∑ i = 1 N i n T – 1 N c e N ŝ o R i, ważonej przez pozostały czas w stosunku do czasu już upłynął, zaokrąglonej do najbliższej liczby całkowitej. Ale jeśli liczba ta była większa niż liczba pacjentów nadal zagrożonych na początku ostatniej przerwy, wybrano tę liczbę. Założenie to formalnie zapisane jest w poniższym równaniu:

wykonujemy krok 2-3.

Krok 7. Następnie wykorzystujemy zgłoszoną całkowitą liczbę zdarzeń, totevents. Obliczamy szacunkową całkowitą liczbę zdarzeń otrzymaną na początku ostatniego przedziału, ∑ k = 1 u p e R n i n t-1 d ^ k . Jeśli jest to większe lub równe totevents Zakładamy, że nie ma więcej zdarzeń lub cenzurowania:

krok 8. Jeśli ∑ k = 1 u p P e R n i n t-1 d ^ k jest mniejsze niż toteventy, ponownie dostosowujemy szacunkową liczbę ocenzurowanych obserwacji w interwale nint, nce N ^ więc R n i N t, przez różnicę w całkowitej liczbie zdarzeń:

następnie ponownie uruchamiamy kroki 2-3, 8 dla ostatniego przedziału, nint, aż szacowana całkowita liczba zdarzeń, ∑ k = 1 u p p e R n i n t – 1 d ^ k, jest równa zgłoszonej całkowitej liczbie zdarzeń, toteventów lub dopóki szacowana całkowita liczba zdarzeń jest mniejsza niż zgłoszona całkowita liczba zdarzeń, ale całkowita liczba cenzurowania w ostatnim przedziale, nce n ^ więc r n i N t, staje się równa zero.

korekta algorytmu dla przypadku „bez numerów zagrożonych”

w tym przypadku istnieje tylko jeden interwał nint = 1. Najpierw Zakładamy, że całkowita liczba ocenzurowana jest równa zero, a następnie postępujemy tak jak w kroku 8.

korekty algorytmu dla przypadku 'brak wszystkich zdarzeń’

w tym przypadku postępujemy tak jak w przypadku 'wszystkie informacje’ z tym wyjątkiem, że nie można dokonać ponownej regulacji przy użyciu całkowitej liczby zdarzeń i dlatego zatrzymujemy się na kroku 6.

dostosowanie algorytmu dla przypadku ” ani ”

gdy nie podano ani całkowitej liczby zdarzeń, ani liczb zagrożonych po rozpoczęciu badania, założyliśmy, że nie było ocenzurowanych obserwacji. Jest to mocne założenie, ale równie silne, jak każde inne założenie, które moglibyśmy przyjąć o cenzurowaniu bez dalszych informacji. Ze względu na brak informacji oczekuje się niższej jakości wyników.

uzyskanie indywidualnych danych pacjenta (IPD) ze zrekonstruowanych danych Kaplana-Meiera

ze zrekonstruowanych parametrów Kaplana-Meiera d ^ k, cê n k, n ^ k dla każdego wyodrębnionego km współrzędnej k = 1,…, N, możemy uzyskać IPD, który wygeneruje te dane. Ten ostatni element kodowania jest w rzeczywistości dość prosty. Za każdym razem, gdy szacowane jest zdarzenie lub cenzura, rejestrowany jest odpowiedni czas, a także wskaźnik zdarzenia (jeden dla zdarzenia i zero Dla cenzury).

ocena odtwarzalności i dokładności

w ćwiczeniu walidacyjnym wykorzystano sześć par krzywych Kaplana-Meiera. Zostały one zaczerpnięte z podzbioru publikacji, które stanowiły część przeglądu look-back metod analizy czasu przeżycia stosowanych w ocenach ekonomicznych . Przeprowadziliśmy rekonstrukcję dwudziestu dwóch prawdopodobieństw przeżycia, siedmiu median czasu przeżycia, sześciu współczynników ryzyka i czterech standardowych błędów log współczynników ryzyka, które zostały zgłoszone w tych czterech publikacjach. Każdy z nich został dwukrotnie zrekonstruowany przez tych samych trzech obserwatorów. Dwóch z trzech obserwatorów nie było zaangażowanych w rozwój algorytmu.

odtwarzalność i dokładność metody oceniono dla każdego z 4 różnych poziomów informacji („wszystkie informacje”, „Brak liczby zagrożonej”, „brak zdarzeń całkowitych” i „żadne”). Aby ocenić różnice między zrekonstruowanymi statystykami a statystykami oryginalnymi, zastosowano skalę naturalną dla prawdopodobieństwa przeżycia, natomiast skalę logarytmiczną dla medianów, HRs i ich niepewności. Krzywe Kaplana Meiera i Cox HRs na podstawie zrekonstruowanych danych oszacowano przy użyciu procedur r survfit i coxph.

wyposażyliśmy standardową dwukierunkową ANOVĘ z powtarzającymi się pomiarami różnic między zrekonstruowanymi wynikami a oryginalnymi wynikami, w skali naturalnej lub logarytmicznej, w zależności od rozważanej statystyki. Składowymi wariancji były przykład, obserwator, przykład × interakcja obserwatora i błąd wewnątrzkomórkowy. Ponieważ wartość P z testu współczynnika F dla interakcji we wszystkich przypadkach przekraczała 10%, połączyliśmy termin interakcji z terminem błędu wewnątrzkomórkowego. Wybrane podejście jest podobne do tego, co określa się w zastosowaniach inżynieryjnych jako „powtarzalność i odtwarzalność miernika”.

odtwarzalność reprezentuje błąd, jeśli pojedynczy obserwator wykonuje pojedynczą rekonstrukcję dla określonej statystyki. Oszacowano to jako sumę błędu wewnątrz-obserwatora i między-obserwatora. W celu uzyskania 95% przedziałów ufności wokół odchyleń standardowych wykorzystano symulację Monte Carlo z zamontowanego modelu ANOVA. Zakłada się, że stopnie swobody dla wewnątrz, pomiędzy i zmiany wyników są zgodne z rozkładami chi-kwadrat. Aby zapewnić solidne wnioskowanie, z każdego z tych rozkładów pobrano 150 000 próbek stopni swobody, tj. dla każdego źródła zmienności. Następnie obliczono średnie szacunki kwadratów, wykorzystując sumę kwadratów uzyskanych przez ANOVA i próbkę uzyskaną przez symulację, dla każdej z 150 000 próbek i dla każdego ze źródeł zmienności. Następnie oszacowano odpowiednie 150 000 odchyleń standardowych w granicach, między nimi i w wyniku, a następnie wyodrębniliśmy percentyle 2,5 i 97,5 w celu uzyskania oszacowań przedziałów ufności.

aby ocenić dokładność, zbadaliśmy średnią różnicę między zrekonstruowanymi statystykami a oryginalnymi. Wynikająca z tego średnia tendencja lub średni błąd (mean error) odzwierciedla systematyczne nadmierne lub niedoszacowanie. 95% przedziały ufności uzyskuje się bezpośrednio z oszacowania odchyleń standardowych podanych przez ANOVA. Zarejestrowaliśmy również błąd absolutny lub średni błąd absolutny (MAE). Ignoruje to kierunek błędów i mierzy ich wielkość, dając miarę absolutnej dokładności zrekonstruowanych wyników. W celu uzyskania 95% przedziałów ufności ponownie zastosowano metodę symulacji, która zakładała, że MEs są normalnie rozłożone. Dla każdej statystyki, aby zapewnić solidne wnioskowanie, pobrano 150 000 próbek z rozkładu normalnego z obserwowaną średnią i wariancją, zgodnie z ANOVA. Następnie obliczyliśmy odpowiednie wartości bezwzględne 150 000 tych liczb i ostatecznie wyodrębniliśmy 2,5 i 97,5 percentyla, aby uzyskać szacunki przedziałów ufności.

wreszcie zanotowaliśmy różnice w różnicy między zrekonstruowanymi i oryginalnymi statystykami, które wynikały z wyboru przykładów, tj. z 22 prawdopodobieństw przeżycia, 7 medianów, 6 godzin i 4 błędów standardowych logów godzin. Daje to dalsze wskazanie na dokładność metody.