Regresja PH Coxa z efektem losowym (frailty)

Kiedy stosować:

Model Coxa z efektem losowym jest rozszerzeniem regresji proporcjonalnych hazardów Coxa na dane, w których obserwacje tworzą grupy: pacjenci leczeni w tym samym ośrodku, zwierzęta z jednego miotu, zdarzenia nawracające u tego samego pacjenta (kolejne infekcje, kolejne hospitalizacje), oczy lub nerki tego samego chorego (zob. wprowadzenie do modeli mieszanych). Czasy przeżycia w obrębie grupy są ze sobą powiązane, bo grupa dzieli nieobserwowane czynniki ryzyka: pacjent ma swoją odporność, ośrodek – swoją jakość opieki. Zwykła regresja Coxa tę zależność pomija, przez co zaniża błędy standardowe współczynników i zawyża istotność.

Efekt losowy, zwany w analizie przeżycia frailty (kruchością), to nieobserwowany czynnik \(w_j\) wspólny dla całej grupy \(j\), który mnoży hazard każdej jej obserwacji:

\[ h_{ij}(t)=w_j\,h_0(t)\,\exp(\beta_1X_{1ij}+\ldots+\beta_kX_{kij}), \]

gdzie:
\(h_{ij}(t)\) – hazard \(i\)-tej obserwacji w grupie \(j\),
\(h_0(t)\) – hazard bazowy,
\(\beta_1,\ldots,\beta_k\) – współczynniki (efekty stałe),
\(w_j\) – frailty grupy \(j\); zakłada się rozkład gamma o średniej 1 i wariancji \(\theta\) (dzielone frailty gamma, ang. shared gamma frailty).

Grupy o \(w_j>1\) mają wyższe ryzyko, niż wynika z ich zmiennych, grupy o \(w_j<1\) – niższe. Wariancja \(\theta\) mówi, jak bardzo grupy różnią się ryzykiem; \(\theta=0\) oznacza brak zależności w grupach i model sprowadza się do zwykłej regresji Coxa. Frailty jest odpowiednikiem losowego wyrazu wolnego z regresji wielorakiej mieszanej, tyle że działa multiplikatywnie na hazard; losowe nachylenia nie są w tym modelu dostępne.

Warunki stosowania

proporcjonalność hazardów, jak w regresji Coxa,
grupy niezależne od siebie; zależność występuje tylko wewnątrz grup,
co najmniej 3 grupy, a do wiarygodnego oszacowania \(\theta\) zwykle kilkanaście lub więcej,
wystarczająca liczba zdarzeń (obserwacji kompletnych); grupy bez zdarzeń wnoszą mało informacji o \(\theta\),
brak silnej współliniowości zmiennych niezależnych.

Uwaga! Dane muszą być przygotowane w układzie długim: każdy wiersz arkusza to jedna obserwacja (czas, status, zmienne niezależne), a osobna kolumna wskazuje grupę (np. identyfikator ośrodka lub pacjenta). Przy zdarzeniach nawracających każdy epizod to osobny wiersz z czasem liczonym od początku tego epizodu.

Definicje miar

Wariancja frailty \(\theta\) – miara zróżnicowania ryzyka między grupami. Frailty ma średnią 1, więc przy \(\theta=0.4\) odchylenie standardowe \(w_j\) wynosi \(\sqrt{0.4}\approx0.63\): typowa grupa ma hazard o 60% wyższy lub niższy, niż wynikałoby z jej zmiennych. Błąd \(\theta\) jest orientacyjny (z profilu wiarygodności) i nie służy do testowania – do tego jest test LR.

\(\tau\) Kendalla – miara zależności czasów przeżycia dwóch obserwacji z tej samej grupy, wynikająca z rozkładu gamma:

\[ \tau=\frac{\theta}{\theta+2}. \]

\(\tau=0\) to niezależność; \(\tau=0.17\) oznacza umiarkowaną zgodność: obserwacja z grupy, w której inna obserwacja przeżyła długo, sama ma większe szanse przeżyć długo.

Test LR: \(\theta=0\) – test ilorazu wiarygodności porównujący model z frailty ze zwykłą regresją Coxa. Hipoteza leży na granicy przestrzeni parametrów, więc wartość \(p\) z rozkładu \(\chi^2_1\) jest dzielona przez dwa (korekta granicy), o czym raport informuje.

Warunkowy hazard ratio \(HR=e^{b}\) – ile razy zmienia się hazard, gdy \(X\) rośnie o jednostkę w obrębie tej samej grupy (przy tym samym frailty). Populacyjny HR ze zwykłej regresji Coxa jest uśredniony po grupach i co do wartości bezwzględnej bliższy 1; nie należy porównywać ich wprost.

Efekty losowe grup – oszacowane \(w_j\) dla każdej grupy oraz \(\log w_j\) z błędem. Grupy z małą liczbą zdarzeń mają \(w_j\) ściągnięte w stronę 1 (nie wiemy o nich dość, by uznać je za nietypowe). Nie testuje się ich istotności; służą do wskazania grup o podwyższonym lub obniżonym ryzyku.

Zintegrowana wiarygodność – logarytm wiarygodności modelu po scałkowaniu po frailty; na jej podstawie liczone są kryteria informacyjne i testy LR (zarówno test \(\theta=0\), jak i test łącznego wpływu zmiennych niezależnych).

Jak ustawić analizę

Statystyka \(\to\) Modele mieszane \(\to\) Regresja PH Coxa z efektem losowym.

Okno ma układ okna regresji Coxa. Wskazujemy Czas przeżycia, Zmienną ucinaną z wartością oznaczającą zdarzenie w polu Kompletne oraz Zmienne niezależne X1, X2, …; zmienne fikcyjne i interakcje przygotowuje się tak samo jak w regresji Coxa. Opcje modelu z efektem losowym:

Zmienna grupująca (ID) – kolumna wskazująca grupę obserwacji (ośrodek, pacjent, miot); liczbowa lub tekstowa. Jest obowiązkowa – bez niej model nie różni się od regresji Coxa,
Efekty losowe grup (frailty) – dołącza do raportu tabelę oszacowanych \(w_j\) dla każdej grupy; przy setkach grup tabela jest długa,
Dołącz reszty i hazard – reszty martyngałowe i dewiancyjne uwzględniające frailty oraz hazard wynikowy \(w_j H_i\) każdej obserwacji; służą do oceny dopasowania i wskazania obserwacji odstających,
Dołącz wykresy – krzywe przeżycia grup (frailty).

Wiązania czasów są rozstrzygane metodą Breslowa, tak jak w regresji Coxa, więc przy \(\theta=0\) wyniki obu analiz są identyczne.

Wyniki i interpretacja

Liczności – liczba obserwacji, kompletnych (zdarzeń) i uciętych oraz liczba grup. Raport przypomina też warunek liczności \(n\ge10k/p\) znany z regresji Coxa.

Dopasowanie modelu – logarytm zintegrowanej wiarygodności modelu, ta sama wielkość dla modelu z samym efektem losowym (bez zmiennych) oraz test ilorazu wiarygodności dla łącznego wpływu zmiennych niezależnych. Kryteria informacyjne AIC, AICc, BIC i pseudo-\(R^2\) jak w regresji Coxa.

Efekt losowy (frailty gamma) – zmienna grupująca, liczba grup, wariancja frailty \(\theta\) z orientacyjnym błędem (profil), \(\tau\) Kendalla oraz test LR: \(\theta=0\) (vs regresja Coxa) z korektą granicy. Istotny wynik oznacza, że struktury grupowej nie wolno pominąć. Podana jest też liczba iteracji EM.

Współczynniki – \(b\), błąd, przedział ufności, statystyka Walda, wartość \(p\) oraz hazard ratio (warunkowy) \(HR=e^{b}\) z przedziałem ufności. Błędy standardowe pochodzą z penalizowanej macierzy informacji i uwzględniają obecność efektu losowego, dlatego są większe niż w zwykłej regresji Coxa.

Funkcje bazowe – skumulowany hazard i przeżycie dla grupy przeciętnej (\(w_j=1\)): przy średnich wartościach zmiennych oraz bazowe, jak w regresji Coxa.

Efekty losowe grup (frailty gamma) – dla każdej grupy: \(w\), \(\log w\) z błędem, liczba zdarzeń, liczność i skumulowany hazard grupy. Grupa z \(w=1.4\) ma hazard o 40% wyższy niż przeciętna grupa o tych samych zmiennych; grupa z \(w=0.6\) – o 40% niższy. Grupy z wieloma zdarzeniami przy małym skumulowanym hazardzie oczekiwanym dostają \(w>1\), grupy bez zdarzeń mimo długiej obserwacji – \(w<1\).

Wykresy – Krzywe przeżycia grup (frailty): krzywa przeżycia dla przeciętnej grupy (\(w=1\)) oraz dla wybranych grup z ich frailty; im wyższe \(w\), tym niżej leży krzywa.

Wykonanie i interpretacja

Przykład 1. (plik nerkiCewnik.pqs)
W badaniu McGilchrista i Aisbett (1991)[190] obserwowano 38 pacjentów dializowanych z założonym cewnikiem. Rejestrowano czas (w dniach) do zakażenia w miejscu wprowadzenia cewnika; po wyleczeniu zakładano cewnik ponownie, więc każdy pacjent ma dwa epizody. Cewnik bywał usuwany z innych przyczyn – takie obserwacje są ucięte. Dane w układzie długim: Pacjent (1–38), Czas (dni), Zakażenie (1 – zakażenie, 0 – obserwacja ucięta), Wiek (lata) i Płeć (0 – mężczyzna, 1 – kobieta). Łącznie 76 obserwacji, 58 zakażeń. Pytamy, czy wiek i płeć wpływają na ryzyko zakażenia i czy pacjenci różnią się indywidualną podatnością.

Dwa epizody tego samego pacjenta nie są niezależne – pacjent ma swoją odporność, higienę, choroby współistniejące. Gdybyśmy policzyli zwykłą regresję Coxa, potraktowalibyśmy 76 epizodów jak 76 różnych osób. Budujemy więc model z efektem losowym: Czas jako czas przeżycia, Zakażenie jako zmienną ucinaną z wartością 1 jako Kompletne, Wiek i Płeć jako zmienne niezależne, Pacjent jako zmienną grupującą. Włączamy Efekty losowe grup (frailty) i Dołącz wykresy.

Krok 1 – czy pacjenci różnią się podatnością. Wariancja frailty \(\theta=0.397\) (błąd 0.234), \(\tau\) Kendalla \(=0.166\). Test LR \(\theta=0\): \(\chi^2=5.21\), \(p=0.011\) (z korektą granicy) – zwykła regresja Coxa nie jest dla tych danych poprawna, pacjenci istotnie różnią się nieobserwowaną podatnością na zakażenie. Odchylenie standardowe frailty \(\sqrt{0.397}\approx0.63\): typowy pacjent ma hazard zakażenia o około 60% wyższy lub niższy, niż wynika z jego wieku i płci.

Krok 2 – dopasowanie. Zintegrowana log-wiarygodność modelu wynosi \(-182.05\), modelu z samym frailty (bez zmiennych) \(-187.95\); test LR łącznego wpływu zmiennych: \(\chi^2=11.78\), \(df=2\), \(p=0.0028\). AIC \(=370.11\).

Krok 3 – współczynniki. Wiek: \(b=0.0055\), \(HR=1.006\) (95% CI: 0.983–1.029), \(p=0.637\) – wiek nie wpływa na ryzyko zakażenia. Płeć: \(b=-1.556\) (błąd 0.445), \(HR=0.211\) (95% CI: 0.088–0.504), \(p=0.0005\). Hazard zakażenia u kobiet jest o 79% niższy niż u mężczyzn – przy porównaniu w obrębie tej samej „podatności”, czyli warunkowo. Błąd współczynnika płci (0.445) jest większy niż dałaby zwykła regresja Coxa, bo model uwzględnia, że informacja z dwóch epizodów tego samego pacjenta częściowo się powtarza.

Krok 4 – pacjenci nietypowi. Tabela efektów losowych grup: pacjent 1 ma \(w=1.44\) (\(\log w=0.36\), błąd 0.48) – dwa zakażenia przy niewielkim oczekiwanym skumulowanym hazardzie (0.63) oznaczają podatność o 44% wyższą od przeciętnej. Pacjent 4 ma \(w=0.57\): również dwa zakażenia, ale przy dużym skumulowanym hazardzie (5.35), czyli po długiej obserwacji – jego podatność jest o 43% niższa od przeciętnej. Błędy \(\log w\) (około 0.5) są duże, bo każdy pacjent ma tylko dwa epizody; frailty pojedynczych pacjentów należy traktować orientacyjnie.

Wykres krzywych przeżycia grup pokazuje krzywą dla przeciętnego pacjenta (\(w=1\)) oraz krzywe pacjentów o wyższym i niższym frailty: pacjent 1 traci „przeżycie bez zakażenia” szybciej, pacjent 4 – wolniej.

\(\theta=0.40\), test LR \(\theta=0\): \(p=0.011\); płeć: \(HR=0.21\) (\(p=0.0005\)); wiek: \(HR=1.01\) (\(p=0.64\))

Pacjenci istotnie różnią się podatnością na zakażenie, a po jej uwzględnieniu kobiety mają hazard zakażenia o 79% niższy niż mężczyźni; wiek nie ma znaczenia.

DLA ZAINTERESOWANYCH

Estymacja: Parametry szacuje się algorytmem EM na profilowanej wiarygodności częściowej (Klein 1992[186], Therneau i Grambsch 2000[194]). W kroku M dopasowywana jest regresja Coxa z przesunięciem (offsetem) \(\log w_j\) i wyznaczany skumulowany hazard bazowy Breslowa \(H_0(t)\); w kroku E frailty grupy jest zastępowane wartością oczekiwaną

\[ w_j=\frac{1/\theta+D_j}{1/\theta+A_j},\qquad A_j=\sum_{i\in j}H_0(t_i)\exp(\beta^TX_i), \]

gdzie \(D_j\) to liczba zdarzeń w grupie, a \(A_j\) jej skumulowany hazard. Wariancja \(\theta\) jest wyznaczana przez maksymalizację zintegrowanej log-wiarygodności

\[ \ell=\sum_{i:\delta_i=1}\left[\eta_i-\log S_0(t_i)\right]+\sum_j\ell_j+D, \]

\[ \ell_j=D_j\log\theta+\log\Gamma(1/\theta+D_j)-\log\Gamma(1/\theta)-(1/\theta+D_j)\log(1+\theta A_j), \]

która przy \(\theta\to0\) przechodzi w wiarygodność częściową Coxa (\(D\) – łączna liczba zdarzeń). Zbieżność przyjmuje się, gdy zmiana \(\ell\) jest mniejsza niż \(10^{-9}\).

Błędy standardowe: z penalizowanej macierzy informacji dla parametrów \((\beta,\log w_1,\ldots,\log w_G)\), gdzie kara wynika z rozkładu gamma; blok \(\beta\) jest odwracany przez dopełnienie Schura.

\(\tau\) Kendalla: dla dzielonego frailty gamma \(\tau=\theta/(\theta+2)\).

Test na granicy: statystyka \(2(\ell_{frailty}-\ell_{Cox})\) ma asymptotycznie rozkład będący mieszanką 50:50 rozkładu zdegenerowanego w zerze i \(\chi^2_1\) (Self i Liang 1987[193]), dlatego wartość \(p\) z \(\chi^2_1\) dzieli się przez dwa.