Tabele przeżycia

Analiza przeżycia pozwala ocenić, jak długo trwa do wystąpienia interesującego badacza zdarzenia – np. jak długo żyją pacjenci po operacji, jak długo trwa remisja choroby, po jakim czasie następuje awaria urządzenia. Poza medycyną (skąd pochodzi nazwa) wykorzystywana jest też w innych dziedzinach, gdzie funkcjonuje pod nazwami: analiza niezawodności (ang. reliability analysis), analiza trafności (ang. duration analysis) czy analiza historii zdarzeń (ang. event history analysis).
Podstawowym narzędziem opisowym są tu tabele trwania życia i krzywe Kaplana-Meiera, pozwalające oszacować, jaka część badanych „przeżywa” do danego momentu czasu. Gdy interesuje nas porównanie czasu przeżycia pomiędzy grupami (np. między różnymi metodami leczenia), stosuje się porównanie krzywych przeżycia. Wpływ wielu zmiennych jednocześnie na czas przeżycia bada się natomiast za pomocą modeli regresji, przede wszystkim regresji proporcjonalnego hazardu Coxa.
Podstawowe pojęcia analizy przeżycia wyjaśnimy na przykładzie długości życia pacjentów po przeszczepie serca:

Zdarzenie \(-\) interesująca badacza zmiana, np. śmierć.
Czas przeżycia \(-\) okres, jaki upłynął między stanem początkowym a wystąpieniem zdarzenia, np. długość życia pacjenta po przeszczepie serca.

Uwaga! Do analizy trzeba przygotować jedną kolumnę z gotowym, wyliczonym czasem. Gdy dysponujemy dwoma datami: początkową i końcową, czas jaki upłynął między nimi wyliczamy przed analizą, korzystając z formuł arkusza.
Obserwacje ucięte (ang. censored) \(-\) obserwacje, dla których nie mamy pełnej informacji o czasie przeżycia. Rozróżniamy tu:
- obserwację kompletną – znana jest zarówno data przeszczepu, jak i data śmierci pacjenta, więc można wyznaczyć dokładny czas przeżycia po przeszczepie;
- obserwację uciętą prawostronnie – nie jest znana data śmierci pacjenta (chory żyje jeszcze w chwili zakończenia badania), więc nie można ustalić dokładnego czasu przeżycia;
- obserwację uciętą lewostronnie – nie jest znana data przeszczepu serca, choć wiadomo, że nastąpił on przed rozpoczęciem badania, więc również w tym przypadku nie można ustalić dokładnego czasu przeżycia.

Uwaga! Zakończenie badania oznacza zakończenie obserwacji danego pacjenta – nie zawsze jest to ten sam moment dla wszystkich badanych. Zakończeniem obserwacji może być np. moment, w którym tracimy kontakt z pacjentem (nie znamy więc jego dalszego czasu przeżycia). Analogicznie, rozpoczęcie badania nie musi dotyczyć wszystkich pacjentów w tym samym czasie.

Tabele przeżycia

Kiedy stosować:

Tabele trwania życia (tabele przeżycia) to najprostszy sposób podsumowania czasu przeżycia badanej grupy – oś czasu dzielimy na równe przedziały (np. co 3 lata) i dla każdego z nich sprawdzamy, ilu badanych do niego dotrwało, ilu zmarło, a ilu zostało utraconych z obserwacji (ucięte). Ich zaletą jest łatwość interpretacji, wadą – konieczność arbitralnego podziału czasu na przedziały.

Definicje miar

Dla każdego przedziału czasu program PQStat wylicza:

liczbę przypadków wchodzących \(-\) liczbę osób, które przeżyły do początku danego przedziału;
liczbę przypadków uciętych \(-\) liczbę osób w danym przedziale zdefiniowanych jako przypadki ucięte;
liczbę przypadków zagrożonych \(-\) liczbę przypadków wchodzących do danego przedziału pomniejszoną o połowę liczby przypadków uciętych w tym przedziale;
liczbę przypadków kompletnych \(-\) liczbę osób, które doświadczyły zdarzenia (np. zmarły) w danym przedziale;
proporcję przypadków kompletnych \(-\) stosunek liczby przypadków kompletnych (zgonów) w danym przedziale do liczby przypadków zagrożonych w tym przedziale;
proporcję przypadków przeżywających \(-\) obliczaną jako 1 minus proporcja przypadków kompletnych w danym przedziale;
skumulowaną proporcję przeżywających (funkcję przeżycia) \(-\) prawdopodobieństwo przeżycia do danego czasu; ponieważ by przeżyć kolejny przedział, trzeba przeżyć wszystkie poprzednie, prawdopodobieństwo to wylicza się jako iloczyn wszystkich poprzednich proporcji przypadków przeżywających

\(\pm\) błąd standardowy funkcji przeżycia;
gęstość prawdopodobieństwa \(-\) oszacowane prawdopodobieństwo doświadczenia zdarzenia (zgonu) w danym przedziale, w przeliczeniu na jednostkę czasu

\(\pm\) błąd standardowy gęstości prawdopodobieństwa;
stopę hazardu \(-\) prawdopodobieństwo (w przeliczeniu na jednostkę czasu), że przypadek, który przeżył do początku danego przedziału, doświadczy zdarzenia (zgonu) właśnie w tym przedziale

\(\pm\) błąd standardowy stopy hazardu.

Uwaga! Gdy w którymś przedziale czasu nie ma żadnej obserwacji kompletnej, program umożliwia zastosowanie korekty – zerowa liczba przypadków kompletnych zastępowana jest wówczas wartością 0.5.

Uzyskane informacje można zilustrować trzema wykresami: wykresem funkcji przeżycia, wykresem gęstości prawdopodobieństwa oraz wykresem stopy hazardu.

Jak ustawić analizę

Okno z ustawieniami opcji tabel trwania życia wywołujemy poprzez menu Statystyka\(\to\)Analiza przeżycia\(\to\)Tabele przeżycia.

Tabele trwania życia tworzy się dla podanych przez badacza przedziałów czasu przeżycia o równej rozpiętości – przedziały te definiuje się, podając ich krok.

Wykonanie i interpretacja

Przykład 1. (plik przeszczep.pqs)
Badano długość życia chorych po przeszczepie wątroby. Przez okres 21 lat poddano obserwacji grupę 89 pacjentów, których wiek w chwili przeszczepu należał do przedziału \(\langle 45 \)lat\(; 60 \)lat\()\). Fragment zebranych danych przedstawia poniższa tabela:

Dane kompletne w tej analizie to te, o których mamy pełną informację o długości życia po przeszczepie, czyli opisane jako „zgon” (są to 53 osoby, co stanowi 59.55%). Dane ucięte to te, o których tej informacji nie mamy, ponieważ w chwili zakończenia badania żyją (jest to 36 osób, co stanowi 40.45%). Budujemy tablice trwania życia tych pacjentów, tworząc 3-letnie przedziały czasowe:

Dla każdego 3-letniego przedziału możemy zinterpretować uzyskane w tabeli wyniki – przykładowo dla osób, które żyją przynajmniej 9 lat od przeszczepu, a zatem weszły do przedziału [9;12):

liczba osób, które przeżyły do 9 lat od przeszczepu to 39;
jest 7 osób, o których wiemy, że żyły przynajmniej 9-12 lat – informacje o nich były zbierane, ale nie wiemy, czy żyły dłużej, gdyż po tym czasie wypadły z badania;
liczba osób zagrożonych zgonem w tym przedziale to 36;
jest 14 osób, o których wiemy, że zmarły od 9 do 12 lat po przeszczepie;
39.4% narażonych pacjentów zmarło 9-12 lat po przeszczepie;
60.6% narażonych pacjentów żyło 9-12 lat po przeszczepie;
procent przeżywających do 9 lat po przeszczepie to 61.4% \(\pm\) 5%;
0.08 \(\pm\) 0.02 – to prawdopodobieństwo zgonu dla każdego roku z przedziału 9-12;
0.16 \(\pm\) 0.04 – to warunkowe, roczne prawdopodobieństwo zgonu dla przypadku, który przeżył do początku przedziału 9-12.

Wyniki przedstawimy na kilku wykresach:

Funkcja przeżycia opada stopniowo, bez gwałtownego załamania

Prawdopodobieństwo przeżycia maleje wraz z upływem czasu, jaki minął od przeszczepu, ale nie obserwujemy okresu, w którym gwałtownie wzrastałoby prawdopodobieństwo zgonu.

Krzywe Kaplana-Meiera

Kiedy stosować:

Krzywe Kaplana-Meiera pozwalają na ocenę czasu przeżycia bez konieczności arbitralnego grupowania obserwacji w przedziały, tak jak dzieje się to w tabelach przeżycia. Estymator ten wprowadzili Kaplan i Meier (1958)[91]. Wykres funkcji przeżycia Kaplana-Meiera ma postać funkcji schodkowej.

Definicje miar

Podobnie jak dla tablic przeżycia, wyliczamy tutaj funkcję przeżycia, czyli prawdopodobieństwo przeżycia do danego czasu. Na podstawie błędu standardowego (formuła Greenwooda) i transformacji logarytmicznej (log-log) budowane są przedziały ufności wokół tej krzywej.
Punkt czasu, przy którym wartość funkcji przyjmuje 0.5, to mediana czasu przeżycia. Mediana ta wskazuje na 50% ryzyko wystąpienia zgonu, czyli przewiduje, że u połowy pacjentów zgon nastąpi w ciągu wskazanego czasu. Zarówno mediana, jak i inne percentyle wyznaczane są jako najkrótszy czas przeżycia, dla którego funkcja przeżycia jest mniejsza lub równa danemu percentylowi. Dla mediany wyznaczany jest przedział ufności metodą „test-based” Brookmeyera i Crowley (1982)[29].
Średnia czasu przeżycia wyznaczana jest jako pole pod krzywą przeżycia. Dane dotyczące czasu przeżycia są zwykle mocno skośne, dlatego w analizie przeżycia mediana jest zwykle lepszą miarą tendencji centralnej niż średnia.

Jak ustawić analizę

Okno z ustawieniami opcji krzywej Kaplana-Meiera wywołujemy poprzez menu Statystyka\(\to\)Analiza przeżycia\(\to\)Analiza Kaplana-Meiera.

Wykonanie i interpretacja

Przykład (1) c.d. (plik przeszczep.pqs)
Przedstawimy długość życia po przeszczepie wątroby przy pomocy krzywej Kaplana-Meiera:

Funkcja przeżycia nie opada gwałtownie zaraz po przeszczepie. Wnioskujemy stąd, że początkowy okres po transplantacji nie jest szczególnie obarczony ryzykiem zgonu.

Mediana czasu przeżycia \(=10\) lat

W ciągu 10 lat od przeszczepu oczekujemy, że zgon nastąpi u połowy pacjentów. Wartość tę zaznaczamy na wykresie, rysując linię w punkcie 0.5 odpowiadającym medianie. W podobny sposób zaznaczamy na wykresie kwartyle.

Przedział ufności dla mediany możemy zobrazować na wykresie, wyrysowując linie pionowe w oparciu o przedział ufności wokół krzywej oraz linię na poziomie 0.5.

Porównywanie krzywych przeżycia

Kiedy stosować:

Funkcje przeżycia można budować oddzielnie dla różnych podgrup, np. oddzielnie dla kobiet i mężczyzn, a następnie porównać je między sobą. Takie porównanie może dotyczyć zarówno dwóch, jak i kilku krzywych jednocześnie, może też badać, czy krzywe układają się w naturalny trend (np. wraz z wiekiem), oraz uwzględniać (korygować) analizę o wpływ dodatkowego czynnika (warstwy).
Porównania \(k\) krzywych przeżycia \(S_1,S_2,\ldots,S_k\) w programie można dokonać przy pomocy trzech testów:

testu log-rank \(-\) najbardziej znanego i szeroko stosowanego, nawiązującego do procedury Mantela-Haenszela dla wielu tabel 2x2 (Mantel-Haenszel 1959[112], Mantel 1966[114], Cox 1972[48]),
uogólnienia Gehana testu Wilcoxona \(-\) wywodzącego się od testu Wilcoxona (Breslow 1970, Gehan 1965[69][70]),
testu Tarone-Ware \(-\) wywodzącego się od testu Wilcoxona (Tarone i Ware 1977[159]).

Wszystkie trzy testy opierają się na tej samej statystyce testowej – różnią je jedynie wagi \(w_j\) przypisywane poszczególnym punktom osi czasu:

test log-rank: \(w_j=1\) – wszystkie punkty osi czasu mają tę samą wagę, co daje większy wpływ na wynik późniejszym wartościom osi czasu;
uogólnienie Gehana testu Wilcoxona: \(w_j=n_j\) – momenty czasowe ważone są liczbą obserwacji w każdym z nich, a więc większą wagę mają wartości początkowe osi czasu;
test Tarone-Ware: \(w_j=\sqrt{n_j}\) – momenty czasowe ważone są pierwiastkiem z liczby obserwacji, co sytuuje ten test pomiędzy dwoma powyższymi.

Warunki stosowania

proporcjonalność hazardu – hazard, definiowany jako nachylenie krzywej przeżycia, jest miarą tego, jak szybko następuje niepożądane zdarzenie; złamanie tego założenia nie dyskwalifikuje całkowicie powyższych testów, ale obniża ich moc, tym bardziej im wcześniej na osi czasu krzywe się przecinają.

Definicje miar

Iloraz Hazardów (ang. Hazard Ratio –\(HR\)) opisuje wielkość obserwowanej różnicy między parą krzywych przeżycia, na podstawie liczby zdarzeń obserwowanych i oczekiwanych w każdej z grup:

\[ HR= \frac{O_1/E_1}{O_2/E_2} \]

Jeśli \(HR>1\) (np. \(HR=2\)), stopień narażenia na niepożądane zdarzenie w pierwszej grupie jest odpowiednio większy niż w grupie drugiej (dwukrotnie dla \(HR=2\)); jeśli \(HR<1\), sytuacja jest odwrotna; \(HR\approx1\) oznacza, że obie grupy są narażone w podobnym stopniu.

Uwaga! Przedział ufności dla \(HR\) wyliczany jest w oparciu o błąd standardowy logarytmu \(HR\) (Armitage i Berry 1994[11]).

Program pozwala też sprawdzić, czy krzywe układają się w trend zgodny z naturalnym porządkiem porównywanych grup (np. kolejnych kategorii wiekowych), oraz porównać krzywe z uwzględnieniem warstw, czyli po skorygowaniu o dodatkowy czynnik mogący wpływać na wynik porównania (np. porównując czas przeżycia względem wieku, oddzielnie w warstwach wyznaczonych przez płeć).

Jak ustawić analizę

Okno z ustawieniami opcji porównania krzywych przeżycia wywołujemy poprzez menu Statystyka\(\to\)Analiza przeżycia\(\to\)Porównanie grup.

Wykonanie i interpretacja

Przykład (1) c.d. (plik przeszczep.pqs)
Różnice dla dwóch krzywych przeżycia
Przeszczepy wątroby dokonywane były w dwóch różnych szpitalach. Sprawdzimy, czy długość życia pacjentów po przeszczepie zależała od szpitala, w którym dokonano przeszczepu. Porównania krzywych przeżycia dla tych szpitali dokonamy w oparciu o wszystkie zaproponowane w programie testy służące temu porównaniu.
Hipotezy:

\[\begin{array}{ll} \mathcal{H}_0: & \text{krzywa przeżycia pacjentów szpitala 1 }=\text{ krzywa przeżycia pacjentów szpitala 2},\\ \mathcal{H}_1: & \text{krzywa przeżycia pacjentów szpitala 1 }\neq \text{ krzywa przeżycia pacjentów szpitala 2}. \end{array}\]

\(p=0.6004\) dla testu log-rank (\(p=0.6959\) dla Gehana, \(p=0.6465\) dla Tarone-Ware); \(HR=1.1499\), 95%CI: \(\langle\)0.6570; 2.0126\(\rangle\)

Brak podstaw do odrzucenia \(\mathcal{H}_0\) – długość życia pacjentów obu szpitali jest podobna (przedział ufności dla \(HR\) zawiera wartość 1).

Różnice dla wielu krzywych przeżycia
Przeszczepy wątroby dokonywane były u ludzi w różnym wieku. Wyróżniono 3 grupy wiekowe: \(\langle 45 \)lat\(; 50 \)lat\()\), \(\langle 50 \)lat\(; 55 \)lat\()\), \(\langle 55 \)lat\(; 60 \)lat\()\). Sprawdzimy, czy długość życia pacjentów po przeszczepie zależy od ich wieku w chwili dokonania przeszczepu.
Hipotezy:

\[\begin{array}{ll} \mathcal{H}_0: & \text{krzywe przeżycia pacjentów w wieku }\langle 45 \text{lat}; 50 \text{lat}), \langle 50 \text{lat}; 55 \text{lat}), \langle 55 \text{lat}; 60 \text{lat})\\ & \text{są podobne,}\\ \mathcal{H}_1: & \text{przynajmniej jedna krzywa przeżycia, spośród powyższych 3 krzywych,}\\ & \text{różni się od pozostałych.}\\ \end{array}\]

\(p=0.0692\) dla testu log-rank (\(p=0.09279\) dla Gehana, \(p=0.0779\) dla Tarone-Ware)

Brak podstaw do odrzucenia \(\mathcal{H}_0\) – długość życia pacjentów z trzech grup wiekowych jest podobna, choć wartości \(p\) są dość bliskie poziomowi istotności 0.05.

Przeglądając wartości hazardu (ilorazu wartości obserwowanych i oczekiwanych niepożądanych zdarzeń), zauważamy, że z każdą kolejną kategorią wiekową są one nieco wyższe: \(\langle\)0.68, 0.93, 1.43\(\rangle\). Choć różnice te nie są istotne statystycznie, mogą sugerować istnienie trendu wzrostu wartości hazardu wraz z wiekiem.
Trend dla kilku krzywych przeżycia
Wprowadzając do testu informację o uporządkowaniu porównywanych kategorii (wykorzystując zmienną wiek, w której przedziały wiekowe ponumerowano odpowiednio 1, 2 i 3), sprawdzimy, czy w porównywanych krzywych istnieje trend.
Hipotezy:

\[\begin{array}{ll} \mathcal{H}_0: & \text{brak trendu w krzywych czasu przeżycia pacjentów po przeszczepie}\\ & \text{(trendu zależnego od wieku pacjentów w chwili przeszczepu),}\\ \mathcal{H}_1: & \text{czym starsi są pacjenci w momencie dokonania przeszczepu, tym większe/mniejsze}\\ & \text{ jest prawdopodobieństwo ich przeżycia określonego odcinka czasu.}\\ \end{array}\]

\(p=0.0237\) dla testu log-rank (\(p=0.0317\) dla Gehana, \(p=0.0241\) dla Tarone-Ware)

Krzywe przeżycia ułożone są w istotnym statystycznie trendzie: czym starszy pacjent w chwili przeszczepu, tym mniejsze prawdopodobieństwo przeżycia określonego odcinka czasu.

Najniżej na wykresie Kaplana-Meiera znajduje się krzywa dla osób w wieku \(\langle\)55 lat; 60 lat), nad nią krzywa dla pacjentów w wieku \(\langle\)50 lat; 55 lat), a najwyżej krzywa dla pacjentów w wieku \(\langle\)45 lat; 50 lat).
Krzywe przeżycia dla warstw
Sprawdzimy teraz, czy obserwowany wcześniej trend jest niezależny od szpitala, w którym dokonano przeszczepu. W tym celu jako zmienną warstwa wybierzemy szpital.

W raporcie najpierw przedstawiona jest analiza poszczególnych warstw, zarówno wyniki testów, jak i wartości hazardu. W warstwie pierwszej widoczny jest trend wzrostu hazardu, choć nieistotny statystycznie; trend o tym samym kierunku (na pograniczu istotności) obserwowany jest w warstwie drugiej. Kumulacja tych trendów we wspólnej analizie warstw pozwoliła uzyskać istotność trendu krzywych przeżycia.

Analiza w warstwach: trend istotny niezależnie od szpitala

Czym starszy pacjent w chwili przeszczepu, tym mniejsze prawdopodobieństwo przeżycia określonego odcinka czasu – niezależnie od szpitala dokonującego przeszczepu.

Analiza porównawcza krzywych przeżycia w korekcji o warstwy daje wynik istotny dla testu log-rank i Tarone-Ware, a nieistotny dla Gehana, co może wskazywać na to, że pojawiające się różnice w krzywych nie są tak widoczne w początkowych okresach czasu przeżycia, jak w okresach późniejszych. Przyglądając się ilorazowi hazardu dla porównywanych parami krzywych

możemy zlokalizować istotne różnice. Najmniejszy iloraz hazardu uzyskano dla porównania krzywej dla najmłodszej grupy z krzywą dla grupy najstarszej: 0.53, 95% przedział ufności \(\langle\)0.26; 1.05\(\rangle\) zawiera wprawdzie wartość 1, ale jest na pograniczu tej wartości, co może sugerować istotne różnice między odpowiadającymi im krzywymi. By potwierdzić to przypuszczenie, można użyć filtru danych w oknie analizy i porównać krzywe parami.

Uwaga! Przy porównaniach parami należy zastosować jedną z poprawek stosowanych przy wielokrotnych porównaniach i zmodyfikować poziom istotności. W tym przypadku dla poprawki Bonferroniego przy trzech porównaniach poziom istotności wyniesie 0.017. Dla uproszczenia posłużymy się tylko testem log-rank.

\(\langle\)45 lat; 50 lat) vs \(\langle\)50 lat; 55 lat)

\(\langle\)45 lat; 50 lat) vs \(\langle\)55 lat; 60 lat)

\(\langle\)50 lat; 55 lat) vs \(\langle\)55 lat; 60 lat)

Zgodnie z oczekiwaniem istotne statystycznie różnice dotyczą tylko krzywych przeżycia dla najmłodszej i najstarszej grupy wiekowej.

DLA ZAINTERESOWANYCH

Różnice w krzywych przeżycia
Hipotezy:

\[\begin{array}{ll} \mathcal{H}_0: & S_1(t)=S_2(t)=...=S_k(t),\text{dla wszystkich }t,\\ \mathcal{H}_1: & \text{nie wszystkie }S_i(t)\text{ są sobie równe}. \end{array}\]

Wzory:

Statystyka testowa (chi-kwadrat):

\[ \chi^2=U'V^{-1}U \]

gdzie:

\(U_i=\sum_{j=1}^{m}w_j(d_{ij}-e_{ij})\),
\(V\) – macierz kowariancji o wymiarach \((k-1)\times(k-1)\):

diagonala: \(\sum_{j=1}^{m}w_j^2\frac{n_{ij}(n_j-n_{ij})d_j(n_j-d_j)}{n^2_j(n_j-1)}\),
poza diagonalą: \(\sum_{j=1}^{m}w_j^2\frac{n_{ij}n_{lj}d_j(n_j-d_j)}{n^2_j(n_j-1)}\),

\(m\) \(-\) liczba momentów czasowych, w których nastąpiło niepożądane zdarzenie (zgon),
\(d_j=\sum_{i=1}^k d_{ij}\) \(-\) obserwowana liczba niepożądanych zdarzeń w \(j\)-tym momencie czasowym,
\(d_{ij}\) \(-\) obserwowana liczba niepożądanych zdarzeń w \(i\)-tej grupie w \(j\)-tym momencie czasowym,
\(e_{ij}=\frac{n_{ij}d_j}{n_j}\) \(-\) oczekiwana liczba niepożądanych zdarzeń w \(i\)-tej grupie w \(j\)-tym momencie czasowym,
\(n_j=\sum_{i=1}^k n_{ij}\) \(-\) liczba narażonych w \(j\)-tym momencie czasowym.

Statystyka ta ma asymptotycznie (dla dużych liczności oczekiwanych) rozkład \(\chi^2\) z \(df=k-1\) stopniami swobody.

Wyznaczoną na podstawie statystyki testowej wartość \(p\) porównujemy z poziomem istotności \(\alpha\): gdy \(p \le \alpha\), odrzucamy \(\mathcal{H}_0\) na rzecz \(\mathcal{H}_1\); gdy \(p > \alpha\), nie ma podstaw, aby odrzucić \(\mathcal{H}_0\).

Trend w krzywych przeżycia
Hipotezy:

\[\begin{array}{ll} \mathcal{H}_0: & \text{w badanej populacji nie istnieje trend w położeniu krzywych }S_1,S_2,...,S_k,\\ \mathcal{H}_1: & \text{w badanej populacji istnieje trend w położeniu krzywych }S_1,S_2,...,S_k. \end{array}\]

Wzory:

Statystyka testowa (chi-kwadrat):

\[ \chi^2=\frac{(c'U)^2}{c'Vc} \]

gdzie \(c=(c_1,c_2,...,c_k)\) – wektor wag dla porównywanych grup informujący o ich naturalnym porządku (najczęściej kolejne liczby naturalne). Statystyka ta ma asymptotycznie (dla dużych liczności oczekiwanych) rozkład \(\chi^2\) z jednym stopniem swobody.

Uwaga! By można było przeprowadzić analizę trendu w krzywych przeżycia, zmienna grupująca musi być zmienną liczbową, w której wartości informują o naturalnym porządku grup – liczby te traktowane są w analizie jako wagi \(c_1,c_2,...,c_k\).

Krzywe przeżycia dla warstw
Hipotezy dla różnic w krzywych przeżycia:

\[\begin{array}{ll} \mathcal{H}_0: & S_1^*(t)=S_2^*(t)=...=S_k^*(t),\text{dla wszystkich }t,\\ \mathcal{H}_1: & \text{nie wszystkie }S_i^*(t)\text{ są sobie równe}. \end{array}\]

Hipotezy dla analizy trendu w krzywych przeżycia:

\[\begin{array}{ll} \mathcal{H}_0: & \text{w badanej populacji nie istnieje trend w położeniu krzywych }S_1^*,S_2^*,...,S_k^*,\\ \mathcal{H}_1: & \text{w badanej populacji istnieje trend w położeniu krzywych }S_1^*,S_2^*,...,S_k^*, \end{array}\]

gdzie \(S_1^*(t), S_2^*(t), ..., S_k^*(t)\) – krzywe przeżycia po korekcji o zmienną wyznaczającą warstwy. Obliczenia dla statystyk testowych bazują na formułach opisanych dla testów nieuwzględniających warstw, z tą różnicą, że macierz \(U\) i \(V\) zastąpiona jest sumą macierzy \(\sum_{l=1}^L U\) i \(\sum_{l=1}^L V\), sumowaną po warstwach utworzonych przez zmienną korygującą (\(l=1,2,...,L\)).

Regresja proporcjonalnego hazardu Coxa

Kiedy stosować:

Regresja Coxa, znana też jako model proporcjonalnego hazardu Coxa (Cox PH model) – Cox D.R. (1972)[48], jest najszerzej stosowanym podejściem regresyjnym w analizie przeżycia. Pozwala zbadać jednoczesny wpływ wielu zmiennych niezależnych (\(X_1\), \(X_2\), \(\ldots\), \(X_k\)) na czas przeżycia, np. wpływ sposobu leczenia, wieku i innych czynników rokowniczych na długość życia. Jest to podejście w pewnym sensie nieparametryczne – nie wymaga znajomości kształtu funkcji hazardu ani przeżycia, co jest źródłem jego popularności.

Warunki stosowania

proporcjonalność hazardu – wpływ każdej zmiennej niezależnej na hazard musi być stały w czasie (zob. Analiza reszt modelu poniżej);
brak silnej współliniowości zmiennych niezależnych;
odpowiednio duża liczba obserwacji względem liczby szacowanych parametrów.

Uwaga! Liczba obserwacji powinna być co najmniej dziesięciokrotnie większa niż stosunek liczby szacowanych parametrów modelu (\(k\)) do mniejszej z proporcji liczności uciętych lub kompletnych (\(p\)), czyli \(n\ge 10k/p\) (Peduzzi i in. 1995[132]).

Uwaga! Zmienne niezależne nie powinny być współliniowe – współliniowość prowadzi do niepewnej estymacji i bardzo wysokich wartości błędów.

Definicje miar

Funkcja, na której oparty jest model, opisuje hazard wynikowy jako iloczyn dwóch wielkości, z których tylko jedna zależy od czasu (\(t\)):

\[ h(t,X_1,X_2,...,X_k)=h_0(t)\cdot \exp\left(\sum_{i=1}^k\beta_iX_i\right), \]

gdzie \(h(t,X_1,...,X_k)\) to wynikowy hazard (zmieniające się w czasie ryzyko, zależne od innych czynników, np. sposobu leczenia), \(h_0(t)\) to hazard bazowy (hazard przy założeniu, że wszystkie zmienne objaśniające są równe zero), a \(\sum_{i=1}^k\beta_iX_i\) to kombinacja (najczęściej liniowa) zmiennych niezależnych \(X_1,...,X_k\) i parametrów modelu \(\beta_1,...,\beta_k\), nazywanych współczynnikami regresji. Omówienie przygotowania zmiennych fikcyjnych i interakcji przedstawiono w rozdziale Przygotowanie zmiennych do analizy w modelach wielowymiarowych Przygotowanie zmiennych do analizy w modelach wielowymiarowych.
Współczynniki szacowane są tzw. „częściową” metodą największej wiarygodności (algorytm iteracyjny Newtona-Raphsona) – „częściową”, ponieważ poszukiwanie maksimum funkcji wiarygodności \(L\) odbywa się wprost tylko dla danych kompletnych, natomiast dane ucięte są uwzględniane pośrednio. Korekcja rang wiązanych oparta jest na metodzie Breslowa[28].

Uwaga! Kryterium zbieżności algorytmu iteracyjnego Newtona-Raphsona można kontrolować dwoma parametrami: limitem liczby iteracji oraz kryterium zbieżności (wartością, poniżej której poprawa estymacji uznawana jest za nieznaczącą, co kończy działanie algorytmu).

Iloraz Hazardu (HR) – dla każdej zmiennej niezależnej wyliczany jest jednostkowy Iloraz Hazardu (ang. Hazard Ratio - \(HR\)):

\[ HR_i=e^{\beta_i}. \]

Wyraża on zmianę ryzyka niepożądanego zdarzenia, gdy dana zmienna niezależna rośnie o 1 jednostkę, przy pozostałych zmiennych modelu utrzymanych na stałym poziomie. Interpretacja:

\(HR >1\) – stymulujący wpływ zmiennej: o ile wzrasta ryzyko wystąpienia zdarzenia, gdy zmienna niezależna rośnie o jeden poziom;
\(HR <1\) – destymulujący wpływ zmiennej: o ile spada ryzyko wystąpienia zdarzenia, gdy zmienna niezależna rośnie o jeden poziom;
\(HR\approx1\) – zmienna nie ma wpływu na wystąpienie zdarzenia.

Uwaga! Jeśli analizę przeprowadzamy dla modelu innego niż liniowy lub uwzględniamy interakcję, wówczas – analogicznie jak w modelu regresji logistycznej – odpowiedni \(HR\) wyliczamy, zmieniając formułę kombinacji zmiennych niezależnych.

Weryfikacja modelu – na podstawie współczynnika oraz jego błędu szacunku wnioskujemy, czy dana zmienna niezależna wywiera istotny wpływ na zmienną zależną; posługujemy się w tym celu testem Walda, weryfikującym hipotezę \(\mathcal{H}_0: \beta_i=0\) (równoważnie \(HR_i=1\)) wobec \(\mathcal{H}_1: \beta_i\ne 0\) (\(HR_i\ne1\)).
Jakość całego zbudowanego modelu ocenia się natomiast kilkoma ogólnymi miarami:

kryteria informacyjne (\(AIC\), \(AICc\), \(BIC\)) – opierają się na entropii informacji niesionej przez model, czyli szacują utraconą informację, gdy dany model opisuje badane zjawisko; wybieramy model o minimalnej wartości danego kryterium. Są rodzajem kompromisu między dobrocią dopasowania a złożonością modelu – rosną wraz z liczbą parametrów modelu, tym szybciej, im mniejsza jest liczba obserwacji. \(AICc\) to poprawka \(AIC\) zalecana przy małych próbach, a \(BIC\) dodatkowo silniej karze złożoność modelu (Volinsky i Raftery 2000[165]). Kryterium informacyjne nie jest miarą absolutną – jeśli wszystkie porównywane modele źle opisują rzeczywistość, w żadnym z kryteriów nie znajdziemy tego ostrzeżenia.
Pseudo R\(^2\) (McFadden R\(^2\)) – miara dopasowania modelu, odpowiednik współczynnika determinacji wielorakiej \(R^2\) z liniowej regresji wielorakiej; przyjmuje wartości z przedziału \(\langle0;1)\), gdzie wartości bliskie 1 oznaczają dobre dopasowanie, a 0 – całkowity brak dopasowania. Ze względu na wrażliwość na liczbę zmiennych w modelu wyznacza się też jego poprawione wersje: \(R^2_{Nagelkerke}\) i \(R^2_{Cox-Snell}\).
test ilorazu wiarygodności – podstawowe narzędzie oceny istotności wszystkich zmiennych modelu łącznie; weryfikuje hipotezę \(\mathcal{H}_0\): wszystkie \(\beta_i=0\) wobec \(\mathcal{H}_1\): istnieje \(\beta_i\ne0\).
AUC zależne od czasu – w analizie przeżycia program wyznacza AUC zależne od czasu (Heagerty i in. 2000[81], Blanche i in. 2013[21]): dla każdego czasu \(t\), w którym wystąpiło zdarzenie badane, obserwacje dzielone są na przypadki (zdarzenie do czasu \(t\)) i kontrole (obserwowane dłużej niż \(t\), bez zdarzenia), a markerem jest kombinacja liniowa zmiennych i parametrów modelu. \(AUC(t)\) to prawdopodobieństwo, że losowo wybrany przypadek ma wyższą wartość markera niż losowo wybrana kontrola; cenzurowanie uwzględniane jest przez ważenie obserwacji odwrotnością prawdopodobieństwa cenzurowania (IPCW). Interpretacja jest taka jak dla klasycznego \(AUC\): 0.5 odpowiada rozróżnianiu losowemu, 1 – doskonałemu, a wartości powyżej ok. 0.7–0.8 uznaje się zwykle za dobrą zdolność dyskryminacyjną. Ponieważ \(AUC(t)\) jest funkcją czasu, w raporcie podawane jest średnie \(AUC(t)\) oraz jego wartości w kwartylach i na końcu czasów zdarzeń, a wykres \(AUC(t)\) względem czasu (opcja Krzywa ROC w oknie Dołącz wykres) pokazuje, czy zdolność modelu do przewidywania zdarzeń jest stabilna dla różnych horyzontów prognozy, czy też zmienia się z czasem. Dla najpóźniejszych czasów, gdy w zbiorze kontroli pozostaje już niewiele obserwacji, estymator \(AUC(t)\) staje się niestabilny – skoków na samym końcu wykresu nie należy nadinterpretować. Dla \(AUC(t)\) nie wyznacza się punktu odcięcia ani testu istotności.

Analiza reszt modelu pozwala zweryfikować założenia modelu – przede wszystkim zlokalizować wartości odstające i ocenić proporcjonalność hazardu (opracowano szereg graficznych metod pozwalających na ocenę adekwatności modelu proporcjonalnego hazardu – Lee i Wang 2003[100]). Ze względu na obserwacje ucięte reszt nie liczy się tu w standardowy sposób (jako różnicę wartości obserwowanej i przewidywanej), lecz jako jeden z trzech rodzajów: Martingale, Deviance lub Schoenfeld, które można wyrysować względem czasu lub zmiennych niezależnych.

Martingale: – interpretowane jako różnica w czasie \([0,t]\) pomiędzy obserwowaną a przewidywaną przez model licznością zdarzeń niepożądanych. Wartość oczekiwana wynosi 0, ale rozkład jest skośny (od \(-\infty\) do 1), co utrudnia interpretację wykresu.
Deviance: – podobnie jak Martingale, asymptotycznie mają wartość oczekiwaną 0, ale przy odpowiednio dopasowanym modelu rozkładają się symetrycznie wokół zera z odchyleniem standardowym 1. Wartość dodatnia oznacza, że badany obiekt przeżył krócej niż przewidywał model, ujemna – że dłużej. To przede wszystkim narzędzie identyfikujące wartości odstające (w raporcie reszty oddalone o ponad 3 odchylenia standardowe od 0 oznaczane są kolorem czerwonym).
Schoenfeld: – wyliczane oddzielnie dla każdej zmiennej niezależnej, tylko dla obserwacji kompletnych; ich suma i wartość oczekiwana wynoszą 0. Pozwalają zidentyfikować zmienną łamiącą założenie proporcjonalności hazardu – jest to ta, której wykres reszt względem czasu układa się w systematyczny wzór (najczęściej bada się liniową zależność reszt od czasu). Równomierne rozłożenie punktów wokół 0 świadczy o spełnieniu założenia proporcjonalności hazardu przez daną zmienną.

Uwaga! Gdy dla którejś zmiennej modelu nie jest spełnione założenie proporcjonalności hazardu, rozwiązaniem może być wykonanie analizy Coxa oddzielnie dla każdego poziomu tej zmiennej (stratyfikacja) lub rozszerzenie modelu o odpowiednią zmienną zależną od czasu.

Jak ustawić analizę

Okno z ustawieniami opcji Regresji Coxa wywołujemy poprzez menu Statystyka\(\to\)Analiza przeżycia\(\to\)Regresja PH Coxa.

DLA ZAINTERESOWANYCH

Błąd standardowy współczynnika

\[ SE_b=\sqrt{diag(H^{-1})_b}, \]

gdzie \(diag(H^{-1})\) to główna przekątna macierzy kowariancji.

Test Walda dla pojedynczego współczynnika
Hipotezy: \(\mathcal{H}_0: \beta_i=0\) wobec \(\mathcal{H}_1: \beta_i\ne 0\) (równoważnie \(\mathcal{H}_0: HR_i=1\) wobec \(\mathcal{H}_1: HR_i\ne 1\)). Wzory: Statystyka testowa:

\[ \chi^2=\left(\frac{b_i}{SE_{b_i}}\right)^2 \]

ma asymptotycznie (dla dużych liczności) rozkład \(\chi^2\) z 1 stopniem swobody.

Kryteria informacyjne (oparte na \(L_{FM}\) – maksimum funkcji wiarygodności modelu pełnego, \(L_0\) – maksimum funkcji wiarygodności modelu z samym wyrazem wolnym, \(k\) – liczbie parametrów, \(d\) – liczbie zdarzeń niepożądanych):

\[ AIC=-2\ln L_{FM}+2k,\qquad AICc=AIC+\frac{2k(k+1)}{d-k-1},\qquad BIC=-2\ln L_{FM}+k\ln(d). \]

Pseudo R\(^2\)

\[ R^2_{Pseudo}=1-\frac{\ln L_{FM}}{\ln L_0},\qquad R^2_{Nagelkerke}=\frac{1-e^{-(2/d)(\ln L_{FM}-\ln L_0)}}{1-e^{(2/d)\ln L_0}}, \qquad R^2_{Cox-Snell}=1-e^{\frac{(-2\ln L_0)-(-2\ln L_{FM})}{d}}. \]

Test ilorazu wiarygodności dla wszystkich zmiennych modelu
Hipotezy: \(\mathcal{H}_0\): wszystkie \(\beta_i=0\) wobec \(\mathcal{H}_1\): istnieje \(\beta_i\neq0\). Wzory: Statystyka testowa:

\[ \chi^2=-2\ln(L_0/L_{FM})=-2\ln(L_0)-(-2\ln(L_{FM})) \]

ma asymptotycznie (dla dużych liczności) rozkład \(\chi^2\) z \(k\) stopniami swobody.

Wyznaczoną na podstawie statystyki testowej wartość \(p\) porównujemy z \(\alpha\): gdy \(p \le \alpha\), odrzucamy \(\mathcal{H}_0\) na rzecz \(\mathcal{H}_1\); gdy \(p > \alpha\), nie ma podstaw, aby odrzucić \(\mathcal{H}_0\).

AUC zależne od czasu
Dla markera \(M=\sum_{i=1}^k b_iX_i\) i czasu \(t\), w którym wystąpiło zdarzenie badane, estymator typu cumulative/dynamic z wagami IPCW (Heagerty i in. 2000[81], Blanche i in. 2013[21]) ma postać:

\[ AUC(t)=\frac{\sum_{i}\sum_{j} w_i(t)\,w_j(t)\left[\mathbf{1}(M_i>M_j)+\frac{1}{2}\mathbf{1}(M_i=M_j)\right]}{\sum_i w_i(t)\sum_j w_j(t)}, \]

gdzie sumowanie po \(i\) obejmuje przypadki (\(T_i\le t\), \(\delta_i=1\)) z wagą \(w_i(t)=1/\hat G(T_i^-)\), a po \(j\) – kontrole (\(T_j>t\)) z wagą \(w_j(t)=1/\hat G(t)\); \(\hat G\) jest estymatorem Kaplana-Meiera rozkładu cenzurowania (za „zdarzenie” uznawane jest w nim ucięcie obserwacji). \(AUC(t)\) wyliczane jest we wszystkich unikalnych czasach zdarzeń badanych, a podawane w raporcie średnie \(AUC(t)\) jest średnią arytmetyczną tych wartości.

Porównywanie modeli regresji PH Coxa

Kiedy stosować:

Ze względu na możliwość jednoczesnej analizy wielu zmiennych niezależnych w jednym modelu regresji Coxa, pojawia się problem wyboru optymalnego modelu. Wybierając zmienne niezależne, należy pamiętać, by w modelu znajdowały się zmienne silnie związane z czasem przeżycia i słabo skorelowane między sobą. Porównywanie modeli z różną liczbą zmiennych pozwala ocenić, czy dodanie (lub usunięcie) danej zmiennej istotnie poprawia model.

Definicje miar

Porównując modele, zwracamy uwagę na kryteria informacyjne (\(AIC\), \(AICc\), \(BIC\)) oraz współczynniki dopasowania modelu (\(R_{Pseudo}^2\), \(R^2_{Nagelkerke}\), \(R^2_{Cox-Snell}\)), a także na wynik testu ilorazu wiarygodności, weryfikującego, czy maksimum funkcji wiarygodności modelu pełnego (\(L_{FM}\)) różni się istotnie od modelu zredukowanego (\(L_{RM}\)).
Decyzję o wyborze modelu podejmujemy na podstawie wszystkich powyższych wielkości łącznie. Jeśli porównywane modele nie różnią się istotnie, powinniśmy wybrać ten z mniejszą liczbą zmiennych – oznacza to bowiem, że zmienne obecne w modelu pełnym, a nieobecne w zredukowanym, nie wnoszą istotnej informacji. Jeśli różnica jest istotna statystycznie, oznacza to, że model z większą liczbą zmiennych (o większym \(R^2\)) jest istotnie lepszy.
W programie PQStat porównywanie modeli można przeprowadzić na dwa sposoby:

ręcznie \(-\) badacz sam buduje model pełny (z większą liczbą zmiennych) i zredukowany (powstały z pełnego po usunięciu zmiennych zbędnych z punktu widzenia badanego zjawiska), a następnie na podstawie wyników porównania wybiera lepszy z nich;
automatycznie \(-\) program sam buduje kolejne modele i wskazuje najlepszy według wybranego kryterium (opis poniżej).

Automatyczny dobór modelu
Składniki modelu. Jednostką, którą metody automatyczne dodają do modelu lub z niego usuwają, jest składnik: pojedyncza zmienna, cała zmienna fikcyjna (wszystkie jej kolumny 0–1 naraz) albo interakcja. Obowiązuje zasada hierarchii: interakcja może znaleźć się w modelu tylko wtedy, gdy są w nim jej efekty główne, i tylko ona może zostać usunięta przed nimi.
Kryterium doboru wybieramy w oknie analizy:

AIC, AICc lub BIC – kryteria informacyjne; lepszy jest model o mniejszej wartości kryterium. Ruch (dodanie lub usunięcie składnika) jest przyjmowany tylko wtedy, gdy obniża kryterium;
\(p\) testu zagnieżdżonego – test ilorazu wiarygodności porównujący model przed zmianą i po zmianie, z zadanym progiem \(p\). Składnik jest usuwany, gdy \(p\) jest większe od progu, a dodawany, gdy \(p\) jest mniejsze od progu. Ta sama wartość \(p\) jest pokazywana w raporcie jako porównanie sąsiednich modeli.

Metody doboru:

eliminacja wstecz – start od modelu ze wszystkimi zmiennymi; w każdym kroku usuwany jest ten składnik, którego usunięcie najbardziej obniża kryterium; koniec, gdy żadne usunięcie już go nie obniża;
selekcja w przód – start od modelu zerowego (bez zmiennych, sam hazard bazowy); w każdym kroku dodawany jest składnik najbardziej obniżający kryterium; koniec, gdy żadne dodanie już go nie obniża;
krokowa dwukierunkowa – start od modelu zerowego; w każdym kroku rozważane są zarówno dodania, jak i usunięcia składników i wybierany jest najlepszy ruch; koniec, gdy żaden ruch nie poprawia kryterium;
wszystkie podzbiory – dopasowywane są wszystkie modele zgodne z zasadą hierarchii (maksymalnie 15 składników), a w raporcie prezentowany jest najlepszy model każdej wielkości; metoda korzysta wyłącznie z kryteriów informacyjnych. Przy wielu zmiennych obliczenia mogą trwać dłużej, bo liczba dopasowywanych modeli rośnie wykładniczo.

Wszystkie modele są liczone na tym samym zbiorze obserwacji – wiersze z brakami w którejkolwiek ze zmiennych kandydujących są odrzucane raz, przed analizą. Dzięki temu kryteria informacyjne i testy ilorazu wiarygodności są między modelami porównywalne.
Raport zawiera: metodę i kryterium doboru, liczbę składników i dopasowanych modeli, kryteria informacyjne modelu zerowego, kolejne modele ścieżki doboru (każdy z wartością \(-2\ln L\), kryteriami \(AIC\), \(AICc\), \(BIC\), współczynnikami \(R_{Pseudo}^2\), \(R^2_{Nagelkerke}\), \(R^2_{Cox-Snell}\) oraz tabelą współczynników i ilorazów hazardu) i porównanie każdej sąsiedniej pary modeli testem ilorazu wiarygodności (statystyka \(\chi^2\), stopnie swobody, wartość \(p\)). Najlepszy model według wybranego kryterium jest oznaczony w raporcie. Jeśli żaden składnik nie poprawia modelu zerowego, raport informuje, że najlepszy jest model zerowy. Przy metodzie wszystkich podzbiorów sąsiednie modele mogą nie być zagnieżdżone – wówczas test nie jest wykonywany.

Uwaga! Wartości \(p\) i przedziały ufności w modelu wybranym automatycznie są zbyt optymistyczne, ponieważ ten sam zbiór danych posłużył do wyboru zmiennych i do oszacowania ich efektów. Dobór automatyczny warto traktować jako pomoc w poszukiwaniu modelu, a ostateczny wybór oprzeć na wiedzy merytorycznej. Wybrany model należy zbudować ponownie w oknie Regresja PH Coxa i zweryfikować założenie proporcjonalności hazardów.

Jak ustawić analizę

Okno z ustawieniami opcji porównywania modeli wywołujemy poprzez menu Statystyka\(\to\)Analiza przeżycia\(\to\)Regresja PH Coxa \(-\) porównywanie modeli.

Wykonanie i interpretacja

Przykład 2. (plik: remisjaBiałaczka.pqs)
Analiza opiera się na danych dotyczących białaczki opisanych w pracy Freireich i innych 1963[65] i analizowanych dalej przez wielu autorów, m.in. Kleinbaum i Klein 2005[94]. Dane zawierają informację o czasie (w tygodniach) pozostawania w remisji, aż do momentu wycofania pacjenta z badania z powodu wyjścia z remisji (nawrotu objawów) lub ucięcia informacji o pacjencie. Wyjście z remisji to zdarzenie niekorzystne – traktowane jest jako obserwacja kompletna. Obserwacja jest ucięta, jeśli pacjent pozostaje w badaniu do końca i remisja nie nastąpi, lub jeśli opuści badanie.
Pacjenci przydzieleni zostali do dwóch grup: leczonej tradycyjnie (oznaczonej jako 1, czasami nazywanej „grupą placebo”) i leczonej nową metodą (oznaczonej jako 0). Zebrano informację o płci pacjentów (1=mężczyzna, 0=kobieta) oraz o wartościach wskaźnika liczby białych krwinek oznaczonego jako „log WBC”, który jest znanym czynnikiem prognostycznym.

Celem badania jest określenie wpływu sposobu leczenia na czas pozostawania w remisji przy uwzględnieniu możliwych czynników wikłających (confounder) i interakcji. W analizie uwagę skupimy na zmiennej „Rx” (1=placebo, 0=nowe leczenie); zmienną „log WBC” umieścimy w modelu jako możliwy czynnik wikłający. By ocenić ewentualny wpływ interakcji „Rx” i „log WBC”, rozważymy także trzecią zmienną będącą ich iloczynem, dołączaną do modelu przyciskiem Interakcje.

Budujemy trzy modele Coxa:

Model A: zawiera tylko zmienną „Rx”
Model B: zawiera zmienną „Rx” i potencjalną zmienną wikłającą „log WBC”
Model C: zawiera zmienną „Rx”, zmienną „log WBC” oraz potencjalny efekt ich interakcji: „Rx \(\times\) log WBC”

Zmienna interakcji „Rx” i „log WBC” zawarta w modelu C jest w nim nieistotna (\(p=0.5103\)) według testu Walda. Podobny wynik uzyskamy, porównując testem ilorazu wiarygodności model C z modelem B w module Regresja Coxa \(-\) porównywanie modeli: wynik nieistotny (\(p=0.5134\)) oznacza, że model C (z interakcją) nie jest istotnie lepszy od modelu B (bez interakcji).

Porównanie modelu C i B: \(p=0.5134\) (nieistotne)

Odrzucamy model C – interakcja „Rx” \(\times\) „log WBC” nie wnosi istotnej informacji. Przechodzimy do porównania modelu B i A.

\(HR\) dla „Rx” w modelu B wynosi 3.65, co oznacza, że hazard dla „grupy placebo” jest około 3.6 razy większy niż dla grupy leczonej nową metodą. Model A, zawierający tylko zmienną „Rx”, nazywany jest zwykle modelem „surowym” (ang. crude model), ponieważ ignoruje efekt potencjalnych zmiennych wikłających – w tym modelu \(HR\) dla „Rx” wynosi już 4.52, czyli sporo więcej niż w modelu B. Ważny jest tu też 95% przedział ufności dla tych estymatorów: dla „Rx” w modelu A ma on szerokość 8.06 (10.09 minus 2.03), a w modelu B jest węższy: 6.74 (8.34 minus 1.60) – model B daje więc bardziej precyzyjną estymację \(HR\) niż model A. By ostatecznie zdecydować, który z nich lepiej oszacowuje efekt leczenia, ponownie porównujemy je testem ilorazu wiarygodności.

Porównanie modelu B i A: \(p<0.0001\) (istotne)

Model B jest istotnie lepszy niż model A – ma najniższe wartości kryteriów informacyjnych (AIC=148.6, AICc=149, BIC=151.4) i wysokie wartości dopasowania (Pseudo \(R^2_{McFadden}=0.2309\), \(R^2_{Nagelkerke}=0.7662\), \(R^2_{Cox-Snell}=0.7647\)).

Dopełnieniem analizy jest przedstawienie dla modelu B krzywych przeżycia (pozostawania w remisji) obu grup leczenia, skorygowanych o wpływ „log WBC”. By je wyrysować, po wybraniu opcji Dołącz wykres zaznaczamy Funkcja przeżycia: w podgrupach..., a następnie – by szybko zbudować wykres dwóch krzywych – wybieramy Szybkie podgrupy i wskazujemy zmienną Rx. Opcja Zaawansowane podgrupy pozwala zbudować dowolną liczbę dowolnie zdefiniowanych krzywych.

Na koniec ocenimy założenia regresji Coxa, analizując reszty modelu względem czasu.

Nie obserwujemy punktów odstających, choć reszty Martingale i Deviance są coraz niższe dla coraz dłuższego czasu. Reszty Schoenfelda mają symetryczny rozkład względem czasu. W przypadku reszt Schoenfelda analizę wykresu można wesprzeć testami sprawdzającymi, czy punkty wykresu reszt układają się w pewien wzór, np. zależność liniową – w tym celu przekopiowujemy reszty Schoenfelda oraz czas do arkusza danych i testujemy szukaną zależność. Wynik takiego testu dla każdej zmiennej wskazuje na spełnienie założenia proporcjonalności hazardu (gdy jest nieistotny statystycznie) lub na jego złamanie (gdy jest istotny) – taką zmienną można wówczas wyłączyć z modelu. Dla zmiennych „log WBC” i „Rx” symetryczny rozkład reszt sugeruje spełnienie założenia proporcjonalności hazardu; potwierdzeniem może być sprawdzenie zależności liniowej Pearsona lub monotonicznej Spearmana dla tych reszt i czasu.
Dalej możemy dodać do modelu zmienną płeć – ostrożnie, bo z literatury wiadomo, że płeć może wpływać na funkcję przeżycia w remisji białaczki w sposób nieproporcjonalny względem czasu. Budujemy więc model Coxa dla trzech zmiennych: „Płeć”, „Rx” i „log WBC”, a przed interpretacją współczynników sprawdzamy reszty Schoenfelda – przekopiowując je wraz z czasem do nowego arkusza i sprawdzając zależność monotoniczną Spearmana. Uzyskujemy \(p=0.0259\) (dla płci), \(p=0.6192\) (dla log WBC) i \(p=0.1490\) (dla Rx), co potwierdza naruszenie założenia proporcjonalności hazardu przez zmienną płeć.

Reszty Schoenfelda vs czas: \(p=0.0259\) dla płci (istotne)

Zmienna „Płeć” łamie założenie proporcjonalności hazardu – model Coxa budujemy oddzielnie dla kobiet i mężczyzn, wykonując analizę dwukrotnie z włączonym filtrem danych (najpierw płeć żeńska, potem męska).

Dla kobiet

Dla mężczyzn

Regresja Coxa ze zmiennymi zależnymi od czasu

Kiedy stosować:

Podstawowym założeniem modelu Coxa jest proporcjonalność hazardu: wpływ każdej zmiennej niezależnej na hazard musi być stały w czasie, czyli Iloraz Hazardów \(HR_i\) nie zmienia się wraz z upływem czasu. Gdy to założenie jest złamane dla którejś ze zmiennych (przesłanką może być np. nielosowy, systematyczny wzór reszt Schoenfelda względem czasu – patrz Analiza reszt modelu), jednym z rozwiązań – obok stratyfikacji analizy (osobne modele budowane dla poszczególnych poziomów tej zmiennej) – jest rozszerzenie modelu o dodatkowy człon opisujący, jak efekt zmiennej zmienia się w czasie.

Definicje miar

Rozszerzenie modelu polega na dodaniu członu będącego iloczynem badanej zmiennej niezależnej \(X_j\) i pewnej funkcji czasu \(g(t)\):

\[ h(t,X_1,...,X_k)=h_0(t)\cdot\exp\left(\sum_{i=1}^k\beta_iX_i+\delta\cdot X_j\cdot g(t)\right), \]

gdzie iloczyn \(X_j\cdot g(t)\) nazywany jest zmienną zależną od czasu (ang. time-dependent covariate). W efekcie hazard cząstkowy związany ze zmienną \(X_j\) przestaje być stały i zależy od czasu:

\[ HR_j(t)=e^{\beta_j+\delta\cdot g(t)}. \]

Test Walda dla współczynnika \(\delta\) pozwala sprawdzić, czy efekt zmiennej \(X_j\) rzeczywiście zmienia się w czasie: \(\delta=0\) oznacza, że efekt \(X_j\) jest stały (założenie proporcjonalności hazardu spełnione), a \(\delta\ne0\) – że efekt zmienia się w czasie (założenie złamane).
W programie PQStat funkcję \(g(t)\) można wybrać spośród czterech dodatkowych pseudo-zmiennych czasu, z którymi buduje się iloczyn z dowolną zmienną niezależną modelu:

czas (liniowo) \(-\) \(g(t)=t\),
ln(czas) \(-\) \(g(t)=\ln(t)\),
sqrt(czas) \(-\) \(g(t)=\sqrt{t}\),
schodkowa(czas) \(-\) \(g(t)\) jest funkcją skokową (Heaviside'a), przyjmującą wartość 0 przed ustalonym progiem czasu i 1 po jego przekroczeniu.

Jak ustawić analizę

Okno z ustawieniami opcji regresji Coxa ze zmiennymi zależnymi od czasu wywołujemy poprzez menu Statystyka\(\to\)Analiza przeżycia\(\to\)Regresja PH Coxa – to samo okno, co dla standardowej regresji Coxa, w którym zaznaczamy dodatkowo opcję Zmienne zależne od czasu. Po jej zaznaczeniu i wskazaniu zmiennej definiującej czas przeżycia, na liście zmiennych do budowy interakcji (przycisk Interakcje) pojawiają się cztery dodatkowe pseudo-zmienne czasu wymienione powyżej – funkcję \(g(t)\) definiujemy więc, wybierając jedną z nich.

Uwaga! Opcja Zmienne zależne od czasu oraz opcja Analiza ryzyk konkurencyjnych (Fine-Gray) wzajemnie się wykluczają.

W oknie ustawień zaawansowanych (przycisk z ikoną koła zębatego), obok znanych już parametrów zbieżności algorytmu Newtona-Raphsona, dostępne są dodatkowo dwa parametry związane ze zmiennymi zależnymi od czasu:

wartość \(a\) używana do wyliczenia \(\ln(t+a)\) dla obserwacji, dla których czas przeżycia wynosi 0 (domyślnie funkcja logarytmiczna nie jest wówczas określona),
próg dla funkcji schodkowej \(-\) może nim być mediana zmiennej definiującej czas (opcja domyślna) albo dowolna inna, wskazana przez badacza wartość.

Wykonanie i interpretacja

Przykład (2) c.d. (plik remisjaBiałaczka.pqs)
W przykładzie dotyczącym remisji białaczki analiza reszt Schoenfelda wykazała naruszenie założenia proporcjonalności hazardu przez zmienną płeć (1=mężczyzna, 0=kobieta; istotna zależność monotoniczna reszt od czasu, p=0.0259), podczas gdy dla zmiennych Rx i log WBC założenie to było spełnione. Zanim zbudujemy model rozszerzony, przyjrzyjmy się, jak to naruszenie wygląda w danych i do czego prowadzi w zwykłym modelu Coxa.
Krok 1. Krzywe Kaplana-Meiera dla płci
Porównanie krzywych przeżycia (pozostawania w remisji) kobiet i mężczyzn w module Porównanie grup daje wynik nieistotny: test log-rank \(\chi^2=0.56\), \(p=0.4555\), iloraz hazardów \(HR_{m/k}=0.78\) (95%CI: \(\langle\)0.38; 1.60\(\rangle\)).

Wykres pokazuje jednak, dlaczego test „nie widzi” różnicy: krzywe przecinają się około 11.–12. tygodnia. W pierwszych tygodniach mężczyźni wychodzą z remisji znacznie szybciej (ich krzywa leży niżej), ale ci, którzy przetrwają ten okres, pozostają w remisji dłużej niż kobiety – po 12. tygodniu krzywa dla mężczyzn leży już wyżej i stabilizuje się na poziomie ok. 0.28, podczas gdy krzywa dla kobiet spada do 0. Przecinające się krzywe to klasyczny obraz złamania założenia proporcjonalności hazardu: hazard mężczyzn względem kobiet nie jest stały, lecz zmienia kierunek w czasie. Test log-rank, który przypisuje wszystkim punktom czasu tę samą wagę, sumuje różnice o przeciwnych znakach – i traci moc.
Krok 2. Model Coxa ignorujący złamanie założenia
Ten sam efekt „znoszenia się” obserwujemy w zwykłym modelu Coxa z trzema zmiennymi: płeć, log WBC i Rx.

Współczynnik przy zmiennej płeć jest nieistotny (\(HR=1.30\), 95%CI: \(\langle\)0.54; 3.14\(\rangle\), \(p=0.5582\)), a wyrysowane z modelu krzywe przeżycia dla obu płci (przy średnich pozostałych zmiennych) leżą blisko siebie i – z założenia modelu – nie mogą się przeciąć. Model uśrednia wczesny, bardzo niekorzystny efekt bycia mężczyzną z późniejszym efektem korzystnym i wychodzi mu „nic”. Wniosek, że płeć nie ma znaczenia dla czasu remisji, byłby więc błędny – jest artefaktem niespełnionego założenia.
Krok 3. Model Coxa ze zmienną zależną od czasu
Zamiast dzielić analizę na dwa oddzielne modele (osobno dla kobiet i dla mężczyzn), budujemy jeden wspólny model, w którym płeć wchodzi w interakcję z funkcją czasu. Ponieważ interesuje nas nie to, czy, ale jak z upływem czasu zmienia się wpływ płci, wybieramy funkcję logarytmiczną \(g(t)=\ln(t)\) – efekt zmiennych klinicznych w analizie przeżycia zmienia się zwykle najszybciej na początku obserwacji i stabilizuje się z czasem, co dobrze oddaje logarytm. W oknie analizy zaznaczamy opcję Zmienne zależne od czasu, a w oknie Interakcje definiujemy iloczyn płeć\(\times\)ln(czas). Model zawiera więc cztery człony: płeć, log WBC, Rx oraz płeć\(\times\)ln(czas).

Współczynnik interakcji płeć\(\times\)ln(czas) wynosi \(\delta=-1.912\) (\(SE=0.922\), \(\chi^2_{Walda}=4.30\), \(p=0.0380\)) i jest istotny statystycznie – potwierdza to formalnie, że efekt płci zmienia się w czasie, czyli to, co wcześniej sugerowały reszty Schoenfelda i przecinające się krzywe Kaplana-Meiera. Istotny stał się też sam współczynnik przy zmiennej płeć: \(\beta_{płeć}=4.057\) (\(p=0.0369\)).
Interpretując te dwa współczynniki, trzeba pamiętać, że \(\beta_{płeć}\) opisuje efekt płci tylko w chwili \(t=1\) (bo \(\ln(1)=0\)), a nie efekt uśredniony w całym okresie obserwacji. Hazard mężczyzn względem kobiet w dowolnym tygodniu \(t\) wyliczamy ze wzoru \(HR_{płeć}(t)=e^{\beta_{płeć}+\delta\cdot\ln(t)}=e^{4.057-1.912\cdot\ln(t)}\):

\(t\) (tygodnie)	1	2	4	8	12	16	23
\(HR_{płeć}(t)\)	57.8	15.4	4.1	1.1	0.50	0.29	0.14

Na początku obserwacji ryzyko wyjścia z remisji jest u mężczyzn wielokrotnie wyższe niż u kobiet, ale przewaga ta szybko maleje; po około 8 tygodniach (\(t^*=e^{4.057/1.912}\approx 8.4\)) hazardy obu płci się zrównują, a później kierunek zależności się odwraca – mężczyźni, którzy pozostali w remisji, mają już ryzyko nawrotu kilkukrotnie niższe niż kobiety. Dokładnie taki obraz pokazywały krzywe Kaplana-Meiera z kroku 1 – model ze zmienną zależną od czasu nie tylko go potwierdza, ale pozwala go zmierzyć, z jednoczesną korekcją o log WBC i Rx.

Model rozszerzony jest też po prostu lepszy od modelu z kroku 2: ma wyższą wiarygodność (\(-2\ln L=138.4\) wobec \(144.2\)), niższe kryteria informacyjne (\(AIC=146.4\) wobec \(150.2\); \(BIC=152.0\) wobec \(154.4\)) i wyższe współczynniki dopasowania (Pseudo \(R^2_{McFadden}=0.264\) wobec \(0.233\)). Test ilorazu wiarygodności porównujący oba modele (różnica \(-2\ln L\) równa \(5.82\) przy 1 stopniu swobody) daje \(p=0.0158\), czyli dodanie członu zależnego od czasu istotnie poprawia model.

\(\delta_{płeć\times\ln(t)}=-1.91\), \(p=0.0380\); \(HR_{płeć}(t)\) od 57.8 (\(t=1\)) przez 1.0 (\(t\approx8.4\)) do 0.14 (\(t=23\))

Efekt płci nie jest stały w czasie: mężczyźni mają znacznie wyższe ryzyko wyjścia z remisji w pierwszych tygodniach, a znacznie niższe później. Zwykły model Coxa (i test log-rank) uśredniają te przeciwne efekty do nieistotnego „braku różnicy”; model ze zmienną zależną od czasu ujawnia i mierzy tę zmianę.

DLA ZAINTERESOWANYCH

Test Walda dla współczynnika interakcji ze zmienną czasu
Hipotezy:

\[\begin{array}{cl} \mathcal{H}_0: & \delta=0 \quad (\text{efekt } X_j \text{ jest stały w czasie \text{-} założenie proporcjonalności hazardu spełnione}),\\ \mathcal{H}_1: & \delta\ne0 \quad (\text{efekt } X_j \text{ zmienia się w czasie \text{-} złamanie założenia proporcjonalności hazardu}). \end{array}\]

Statystyka testowa oraz sposób wnioskowania są analogiczne jak dla testu Walda pozostałych współczynników modelu Coxa.

Analiza ryzyk konkurencyjnych – model Fine'a-Graya

Kiedy stosować:

W wielu badaniach obserwowane jednostki mogą doświadczyć jednego z kilku, wzajemnie wykluczających się rodzajów zdarzeń niepożądanych. Przykładowo pacjent onkologiczny może umrzeć z powodu choroby podstawowej, ale może też umrzeć z zupełnie innej przyczyny (np. sercowo-naczyniowej), zanim dojdzie do zdarzenia będącego przedmiotem badania. Wystąpienie takiego alternatywnego zdarzenia nazywane jest zdarzeniem konkurencyjnym (ang. competing risk) – uniemożliwia ono zaobserwowanie zdarzenia badanego, a jednocześnie nie jest zwykłym cenzurowaniem, ponieważ dostarcza dodatkowej informacji (wiemy, dlaczego obserwacja się zakończyła).
Zastosowanie standardowego modelu Coxa, w którym zdarzenia konkurencyjne potraktowane zostałyby jak zwykłe cenzurowanie, prowadzi do obciążonych oszacowań – standardowa funkcja \(1-S(t)\) przecenia wówczas prawdziwe ryzyko zdarzenia badanego, ponieważ część osób nigdy by go nie doświadczyła (zmarły wcześniej z innej przyczyny). Właściwym narzędziem analizy w takiej sytuacji jest skumulowana funkcja zapadalności (ang. cumulative incidence function, CIF) oraz model regresji Fine'a-Graya (Fine i Gray 1999)[63], oparty na tzw. subdystrybucyjnym hazardzie (ang. subdistribution hazard).

Definicje miar

Skumulowana funkcja zapadalności (CIF) wyraża, dla zdarzenia badanego, prawdopodobieństwo jego wystąpienia do czasu \(t\), przy uwzględnieniu faktu, że zdarzenia konkurencyjne mogą temu zapobiec. Suma CIF dla zdarzenia badanego, CIF dla zdarzeń konkurencyjnych oraz funkcji przeżycia \(S(t)\) w każdym punkcie czasu wynosi 1.
Model subdystrybucyjnego hazardu Fine'a-Graya definiuje subdystrybucyjny hazard \(h^*(t)\) analogicznie do tego, jak model Coxa opisuje standardowy hazard, opiera się jednak na innym niż standardowo zbiorze ryzyka: obiekty, u których wystąpiło zdarzenie konkurencyjne, nie są usuwane ze zbioru ryzyka – pozostają w nim (z odpowiednią wagą) aż do końca obserwacji, symulując sytuację, w której nigdy nie mogłyby doświadczyć zdarzenia badanego, ale nadal są obserwowane. Wagi te (ang. inverse probability of censoring weighting, IPCW) wyliczane są w oparciu o odwrotność estymatora Kaplana-Meiera rozkładu cenzurowania. Parametry modelu szacowane są – podobnie jak w standardowym modelu Coxa – metodą największej wiarygodności (algorytm iteracyjny Newtona-Raphsona), z tą różnicą, że funkcja wiarygodności uwzględnia opisane powyżej wagi IPCW.
Subdystrybucyjny Iloraz Hazardów (sHR) – dla każdej zmiennej niezależnej wyliczany jest subdystrybucyjny iloraz hazardów (ang. subdistribution Hazard Ratio –\(sHR\)). W odróżnieniu od standardowego \(HR\) modelu Coxa (opisującego wpływ zmiennej na tempo występowania zdarzenia badanego wyłącznie u osób nadal „zagrożonych” tym zdarzeniem), \(sHR\) opisuje wpływ zmiennej niezależnej bezpośrednio na skumulowaną funkcję zapadalności \(CIF(t)\), co czyni go łatwiejszym do interpretacji klinicznej (np. do prognozowania) w kontekście ryzyk konkurencyjnych. Kierunek interpretacji (\(sHR>1\), \(sHR<1\), \(sHR\approx1\)) jest analogiczny jak dla standardowego \(HR\) opisanego wcześniej, a istotność poszczególnych współczynników weryfikowana jest również testem Walda.

Uwaga! Model Fine'a-Graya jest osobnym modelem od standardowego modelu Coxa liczonego na tych samych danych, w którym zdarzenia konkurencyjne po prostu wykluczono lub pominięto (uznano za brak zdarzenia) – oba modele odpowiadają na różne pytania badawcze i mogą prowadzić do odmiennych, a czasem wręcz przeciwstawnych wniosków (opisywany w literaturze efekt zmiany kierunku, ang. sign reversal) dla tej samej zmiennej niezależnej.

Jak ustawić analizę

Okno z ustawieniami opcji modelu Fine'a-Graya wywołujemy poprzez menu Statystyka\(\to\)Analiza przeżycia\(\to\)Regresja PH Coxa – to samo okno, co dla standardowej regresji Coxa, w którym zaznaczamy dodatkowo opcję Analiza ryzyk konkurencyjnych (Fine-Gray). Po jej zaznaczeniu, na podstawie wartości zmiennej wskazanej jako zmienna ucinana, pojawia się lista Zdarzenie(a) konkurencyjne:, na której wskazujemy – spośród wartości innych niż kod zdarzenia badanego – jedną lub więcej wartości oznaczających zdarzenia konkurencyjne. Wartości statusu niewybrane ani jako zdarzenie badane, ani jako zdarzenie konkurencyjne pozostają prawdziwym cenzurowaniem.

Uwaga! Z uwagi na inny sposób liczenia reszt modelu (reszty Martingale, Deviance i Schoenfelda nie są zdefiniowane dla CIF), opcja Dołącz reszty i hazard jest dla modelu Fine'a-Graya niedostępna, podobnie jak jednoczesne zaznaczenie opcji Zmienne zależne od czasu. Dostępna jest natomiast opcja Krzywa ROC, czyli AUC zależne od czasu – w wersji uwzględniającej zdarzenia konkurencyjne (patrz niżej).

W raporcie z analizy Fine'a-Graya liczności obserwacji kompletnych i uciętych są rozbite dodatkowo o informację, jaki kod (wartość zmiennej statusu) im odpowiadał, a liczba obserwacji uciętych jest dodatkowo podzielona na poszczególne zdarzenia konkurencyjne. Kolumna standardowego \(HR\) jest zaś zastąpiona kolumną \(sHR\) (Fine-Gray).
Analogicznie do standardowego modelu Coxa (opcja Dołącz wykres), dla modelu Fine'a-Graya program pozwala wyrysować wykresy skumulowanej funkcji zapadalności (CIF):

bazową CIF \(-\) skumulowaną zapadalność przy średnich wartościach kowariantów,
CIF dla wybranych wartości zmiennych \(-\) analogicznie jak funkcję przeżycia w standardowym modelu Coxa, dla wskazanych przez badacza wartości zmiennych niezależnych,
CIF dla podgrup \(-\) oddzielne krzywe CIF budowane dla poszczególnych kategorii wskazanej zmiennej grupującej (opcje Szybkie podgrupy / Zaawansowane podgrupy).

Uwaga! W przeciwieństwie do funkcji przeżycia Kaplana-Meiera, krzywa CIF nie dąży zwykle do 0 ani do 1 – jej wartość graniczna to szacowane prawdopodobieństwo, że dana jednostka kiedykolwiek doświadczy zdarzenia badanego, zanim dojdzie u niej do zdarzenia konkurencyjnego.

AUC zależne od czasu dla ryzyk konkurencyjnych. Zdolność modelu Fine'a-Graya do rozróżniania pacjentów, u których wystąpi zdarzenie badane, oceniamy – tak jak dla standardowego modelu Coxa – za pomocą AUC zależnego od czasu, z jedną istotną różnicą w definicji kontroli. Przypadkami w chwili \(t\) są wyłącznie osoby, które do czasu \(t\) doświadczyły zdarzenia badanego; kontrolami są natomiast zarówno osoby obserwowane dłużej niż \(t\) bez żadnego zdarzenia, jak i osoby, u których do czasu \(t\) wystąpiło zdarzenie konkurencyjne – one bowiem z pewnością nie doświadczą już zdarzenia badanego, a dobry model powinien przypisywać im niższe ryzyko. Interpretacja \(AUC(t)\) pozostaje bez zmian: 0.5 – rozróżnianie losowe, 1 – doskonałe; wykres \(AUC(t)\) względem czasu (opcja Krzywa ROC w oknie Dołącz wykres) pokazuje, jak zmienia się jakość prognozy zdarzenia badanego dla kolejnych horyzontów czasu.

Wykonanie i interpretacja

Przykład 3. (plik: szpiczak_konkurencyjne.pqs)
W badaniu obejmującym 445 pacjentów ze zdiagnozowanym szpiczakiem mnogim porównywano skuteczność nowego leku (zmienna int: 0=leczenie standardowe, 1=nowy lek) w wydłużaniu czasu przeżycia (zmienna czas, w dniach). U pacjentów odnotowywano również wyjściową liczbę leukocytów (wbc) oraz obecność niewydolności nerek (nerki) – powikłania często towarzyszącego szpiczakowi mnogiemu i uznawanego za niekorzystny czynnik rokowniczy.

Zgon pacjenta mógł nastąpić z powodu progresji choroby nowotworowej (zdarzenie badane) albo z zupełnie innej przyczyny, np. sercowo-naczyniowej lub infekcyjnej, niezwiązanej bezpośrednio ze szpiczakiem (zdarzenie konkurencyjne). Informację tę koduje zmienna smierc_konk: 0 – pacjent żyje na koniec obserwacji (cenzurowanie), 1 – zgon z powodu szpiczaka, 2 – zgon z innej przyczyny. W danych dostępna jest też uproszczona zmienna smierc (zgon z jakiejkolwiek przyczyny, bez rozróżnienia powodu), wykorzystywana zwykle do standardowego modelu Coxa.
Model standardowy (nieprawidłowy dla tych danych)
Dla porównania zbudujemy najpierw standardowy model Coxa oparty na zmiennej smierc, w którym zgon z innej przyczyny niż szpiczak (zdarzenie konkurencyjne) błędnie potraktowany jest tak samo jak zgon z powodu szpiczaka (zdarzenie badane):

Model Fine'a-Graya
Następnie budujemy właściwy dla tych danych model Fine'a-Graya, wskazując jako zmienną ucinaną smierc_konk, jako zdarzenie badane kod 1 (zgon z powodu szpiczaka) oraz zaznaczając na liście Zdarzenie(a) konkurencyjne: kod 2 (zgon z innej przyczyny).

Uzyskany \(sHR\) dla zmiennej int wynosi \(sHR=0.67\) (95%CI: \(\langle\)0.25; 1.80\(\rangle\)), \(p=0.4235\), co oznacza brak podstaw, by twierdzić, że nowy lek istotnie wpływa na ryzyko zgonu z powodu szpiczaka mnogiego, gdy poprawnie uwzględni się konkurencyjne ryzyko zgonu z innej przyczyny. Inny obraz daje model standardowy (błędnie ignorujący rozróżnienie przyczyn zgonu), w którym \(HR\) dla int wynosi \(HR=2.18\), \(p=0.008034\) – sugerując istotny, niekorzystny wpływ nowego leku na ryzyko zgonu z jakiejkolwiek przyczyny. To dobra ilustracja opisanego wcześniej efektu zmiany kierunku (sign reversal): lek być może nieznacznie zwiększa ryzyko zgonu z innych przyczyn (np. w wyniku działań niepożądanych), co w modelu standardowym miesza się z ryzykiem zgonu z powodu samej choroby nowotworowej, sugerując pozorny, niekorzystny wpływ leku na przebieg szpiczaka – podczas gdy model Fine'a-Graya pokazuje, że lek nie ma istotnego wpływu na ryzyko zgonu specyficznie z powodu szpiczaka.

\(HR_{int}=2.18\), \(p=0.008034\) (model standardowy) wobec \(sHR_{int}=0.67\), \(p=0.4235\) (model Fine'a-Graya)

Nowy lek pozornie zwiększa ogólne ryzyko zgonu, ale nie zwiększa istotnie ryzyka zgonu specyficznie z powodu szpiczaka – efekt widoczny w modelu standardowym znika po poprawnym uwzględnieniu zdarzeń konkurencyjnych.

Odwrotną sytuację obserwujemy dla niewydolności nerek (nerki) – tu oba modele zgodnie wskazują na silny, istotny statystycznie związek z ryzykiem zgonu, przy czym w modelu Fine'a-Graya efekt jest jeszcze wyraźniejszy: \(sHR=38.65\) (95%CI: \(\langle\)19.47; 76.72\(\rangle\)), \(p<0.000001\), wobec \(HR=11.39\) (95%CI: \(\langle\)7.08; 18.32\(\rangle\)), \(p<0.000001\) w modelu standardowym. Niewydolność nerek jest więc nie tylko istotnym czynnikiem ryzyka zgonu z jakiejkolwiek przyczyny, ale przede wszystkim silnie – i jeszcze silniej, niż sugerowałby model standardowy – związana z ryzykiem zgonu z powodu progresji samego szpiczaka, co jest zgodne z jej dobrze znaną rolą niekorzystnego czynnika rokowniczego w tej chorobie (uszkodzenie nerek towarzyszy zwykle zaawansowanej postaci szpiczaka).

\(sHR_{nerki}=38.65\), 95%CI: \(\langle\)19.47; 76.72\(\rangle\), \(p<0.000001\) (model Fine'a-Graya) wobec \(HR_{nerki}=11.39\), \(p<0.000001\) (model standardowy)

Niewydolność nerek jest silnie i konsekwentnie związana z ryzykiem zgonu z powodu szpiczaka – efekt jest jeszcze wyraźniejszy po poprawnym uwzględnieniu konkurencyjnego ryzyka zgonu z innej przyczyny.

Na zakończenie porównamy krzywe CIF zgonu z powodu szpiczaka w grupach wyróżnionych przez zmienną nerki:

Krzywa CIF dla pacjentów z niewydolnością nerek (nerki=1) rośnie znacznie szybciej i osiąga znacznie wyższe wartości niż krzywa dla pacjentów bez tego powikłania (nerki=0) – już po około 100 dniach obserwacji skumulowane prawdopodobieństwo zgonu z powodu szpiczaka przekracza w tej grupie 0.4, a pod koniec obserwacji zbliża się do 1. U pacjentów bez niewydolności nerek skumulowane prawdopodobieństwo zgonu z powodu szpiczaka pozostaje przez cały okres obserwacji niskie (poniżej 0.15). Wykres potwierdza więc liczbowy wynik \(sHR\): niewydolność nerek wiąże się w tych danych z istotnie wyższym ryzykiem zgonu z powodu szpiczaka mnogiego.
Na koniec ocenimy zdolność zbudowanego modelu Fine'a-Graya do rozróżniania pacjentów, którzy umrą z powodu szpiczaka, wybierając w oknie Dołącz wykres opcję Krzywa ROC, czyli wykres AUC zależnego od czasu:

Wartości \(AUC(t)\) są przez cały okres obserwacji wysokie. W pierwszych tygodniach \(AUC(t)\) rośnie od około 0.91 do około 0.96 i utrzymuje się na tym poziomie przez pierwsze ok. 250 dni, a następnie – w miarę jak wśród kontroli ubywa pacjentów bez zdarzenia i przybywa zgonów z innych przyczyn – obniża się do poziomu 0.88–0.92, na którym pozostaje niemal do końca obserwacji. Oznacza to, że model (a w praktyce przede wszystkim zmienna nerki) bardzo dobrze odróżnia pacjentów, którzy umrą z powodu szpiczaka do danego czasu \(t\), od tych, którzy do tego czasu z tego powodu nie umrą – i to niezależnie od horyzontu prognozy. Gwałtowny spadek \(AUC(t)\) w ostatnim punkcie (ok. 1080 dni) jest artefaktem estymacji: w tak późnym czasie w zbiorze ryzyka pozostaje już tylko kilka obserwacji, a estymator oparty na garstce par przypadek–kontrola staje się niestabilny; wartości z samego końca osi czasu nie należy więc interpretować.

\(AUC(t)\approx 0.88\)–\(0.96\) dla \(t\) od kilku do ok. 1000 dni

Model Fine'a-Graya bardzo dobrze i stabilnie w czasie rozróżnia pacjentów zagrożonych zgonem z powodu szpiczaka – jakość prognozy nie zależy istotnie od przyjętego horyzontu czasu.

DLA ZAINTERESOWANYCH

Skumulowana funkcja zapadalności (CIF)

\[ CIF(t)=P(T\le t,\ \varepsilon=1)=\sum_{t_i\le t}S(t_{i-1})\cdot h_1(t_i), \]

gdzie \(S(t_{i-1})\) to ogólna (uwzględniająca wszystkie przyczyny łącznie) funkcja przeżycia tuż przed \(t_i\), a \(h_1(t_i)\) to hazard specyficzny dla zdarzenia badanego w \(t_i\).

Model subdystrybucyjnego hazardu Fine'a-Graya

\[ h^*(t,X_1,...,X_k)=h^*_0(t)\cdot\exp\left(\sum_{i=1}^k\beta_iX_i\right) \]

Wagi IPCW wykorzystywane do utrzymania w zbiorze ryzyka obiektów po wystąpieniu zdarzenia konkurencyjnego wyliczane są w oparciu o odwrotność estymatora Kaplana-Meiera rozkładu cenzurowania \(\hat G(t)\):

\[ w_i(t)=\frac{\hat G(t)}{\hat G(\min(T_i,t))}, \]

gdzie estymator \(\hat G\) liczony jest tak jak estymator Kaplana-Meiera funkcji przeżycia, ale z zamienionymi rolami – za „zdarzenie” uznawane jest tu prawdziwe cenzurowanie, natomiast zarówno zdarzenie badane, jak i zdarzenia konkurencyjne traktowane są jako cenzurowanie tego pomocniczego estymatora.

Subdystrybucyjny Iloraz Hazardów (sHR)

\[ sHR_i=e^{\beta_i}. \]

AUC zależne od czasu dla ryzyk konkurencyjnych
Formuła jest taka sama jak dla standardowego modelu Coxa (patrz Regresja proporcjonalnego hazardu Coxa), zmienia się jedynie zbiór kontroli i ich wagi: obok osób z \(T_j>t\) (waga \(1/\hat G(t)\)) kontrolami są także osoby, u których do czasu \(t\) wystąpiło zdarzenie konkurencyjne (\(T_j\le t\), \(\varepsilon_j=2\)), z wagą \(1/\hat G(T_j^-)\). Estymator \(\hat G\) jest tym samym estymatorem Kaplana-Meiera rozkładu cenzurowania, który służy do wyznaczania wag IPCW modelu (Blanche i in. 2013[21]).