PQStat - Baza Wiedzy

Spis treści

Regresja proporcjonalnego hazardu Cox'a

Regresja proporcjonalnego hazardu Cox'a

Okno z ustawieniami opcji Regresji Cox'a wywołujemy poprzez menu Statystyka→Analiza przeżycia→Regresja PH Cox'a

Regresja Cox'a, znana również jako model proporcjonalnego hazardu Cox'a (Cox PH model) - Cox D.R. (1972)¹⁾, jest najszerzej stosowanym podejściem regresyjnym w analizie przeżycia. Pozwala na zbadanie wpływu wielu zmiennych niezależnych ( $X_1$ , $X_2$ , $\ldots$ , $X_k$ ) na czas przeżycia. Jest to podejście w pewnym sensie nieparametryczne, a więc obarczone niewieloma założeniami - stąd też wynika jego popularność. Nie musi być znana natura ani kształt funkcji hazardu czy przeżycia, a najważniejszy warunek to założenie, które dotyczy również większości parametrycznych modeli przeżycia czyli proporcjonalność hazardu.

Funkcja, na której oparty jest model proporcjonalnego hazardu Cox'a opisuje hazard wynikowy i jest produktem dwóch wielkości, z których tylko jedna jest zależna od czasu ( $t$ ):

$\begin{displaymath} h(t,X_1,X_2,...,X_k)=h_0(t)\cdot \exp\left(\sum_{i=1}^k\beta_iX_i\right), \end{displaymath}$ gdzie:

[ $h(t,X_1,X_2,...,X_k)$ ] - wynikowy hazard opisujący zmieniające się w czasie ryzyko zależne od innych czynników np. sposobu leczenia,
[ $h_0(t)$ ] - hazard bazowy, czyli hazard przy założeniu, że wszystkie zmienne objaśniające są równe zero,

[ $\sum_{i=1}^k\beta_iX_i$ ] - kombinacja (najczęściej liniowa) zmiennych niezależnych i parametrów modelu,

$X_1,X_2,\ldots X_k$ - zmienne objaśniające, niezależne od czasu,

$\beta_1,\beta_2,\ldots \beta_k$ - parametry.

Zmienne fikcyjne i interakcje w modelu

Omówienie przygotowania zmiennych fikcyjnych i interakcji przedstawiono w rozdziale Przygotowanie zmiennych do analizy w modelach wielowymiarowych.

Korekcja rang wiązanych w regresji Cox'a oparta jest na metodzie Breslow'a²⁾

Model można sprowadzić do postaci liniowej:

$\begin{displaymath} \ln\left(\frac{h(t,X_1,X_2,...,X_k)}{h_0(t)}\right)=\sum_{i=1}^k\beta_iX_i. \end{displaymath}$

Rozwiązaniem równania jest wówczas wektor ocen parametrów $\beta_1,\ldots,\beta_k$ nazywanych współczynnikami regresji:

$b=\left( \begin{array}{ccc} b_1\\ b_2\\ \vdots\\ b_k \end{array} \right).$

Współczynniki te szacowane są poprzez tzw. „częściową” metodę największej wiarygodności. Jest to metoda „częściowa” ponieważ poszukiwanie maksimum funkcji wiarygodności $L$ (w programie użyto algorytm iteracyjny Newton-Raphson) odbywa się tylko dla danych kompletnych, dane ucięte są w tym algorytmie uwzględniane, ale nie bezpośrednio.

Każdy współczynnik obarczony jest pewnym błędem szacunku. Wielkość tego błędu wyliczana jest ze wzoru:

$\begin{displaymath} SE_b=\sqrt{diag(H^{-1})_b} \end{displaymath}$

gdzie:

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

Uwaga! Budując model należy pamiętać, że liczba obserwacji powinna być dziesięciokrotnie większa lub równa stosunkowi liczby szacowanych parametrów modelu ( $k$ ) i mniejszej z proporcji liczności uciętych lub kompletnych ( $p$ ), czyli ( $n\ge 10k/p$ ), Peduzzi P. i inni (1995)³⁾.

Uwaga! Budując model należy pamiętać, że zmienne niezależne nie powinny być współliniowe. W przypadku gdy występuje współliniowość, estymacja może być niepewna a uzyskane wartości błędów bardzo wysokie.

Uwaga! Kryterium zbieżności funkcji algorytmu iteracyjnego Newtona-Raphsona można kontrolować przy pomocy dwóch parametrów: limitu iteracji zbieżności (podaje maksymalną ilość iteracji w jakiej algorytm powinien osiągnąć zbieżność) i kryterium zbieżności (podaje wartość poniżej której uzyskana poprawa estymacji uznana będzie za nieznaczną i algorytm zakończy działanie).

Przykład c.d. (plik: remisjaBiałaczka.pqs)

Iloraz Hazardu (HR)

Dla każdej zmiennej niezależnej wyliczany jest jednostkowy Iloraz Hazardu (ang. Hazard Ratio - $HR$ ): $\begin{displaymath} HR_i=e^{\beta_i}. \end{displaymath}$ Wyraża on zmianę ryzyka niepożądanego zdarzenia, gdy zmienna niezależna rośnie o 1 jednostkę. Wynik ten jest skorygowany o pozostałe zmienne niezależne znajdujące się w modelu w ten sposób, że zakłada iż pozostają one na stałym poziomie podczas, gdy badana zmienna niezależna rośnie o jednostkę.

Wartość $HR$ interpretujemy następująco:

$HR >1$ oznacza stymulujący wpływ badanej zmiennej niezależnej na wystąpienie niepożądanego zdarzenia, tj. mówi o ile wzrasta ryzyko na wystąpienie niepożądanego zdarzenia, gdy zmienna niezależna wzrasta o jeden poziom.
$HR <1$ oznacza destymulujący wpływ badanej zmiennej niezależnej na wystąpienie niepożądanego zdarzenia, tj. mówi o ile spada ryzyko na wystąpienie niepożądanego zdarzenia, gdy zmienna niezależna wzrasta o jeden poziom.
$HR\approx1$ oznacza, że badana zmienna niezależna nie ma wpływu na wystąpienie niepożądanego zdarzenia.

Uwaga!

Jeśli analizę przeprowadzamy dla modelu innego niż liniowy, lub uwzględniamy interakcję, wówczas analogicznie jak w modelu regresji logistycznej na podstawie ogólnego wzoru możemy wyliczyć odpowiedni $HR$ zmieniając formułę będącą kombinacją zmiennych niezależnych.

Przykład c.d. (plik: remisjaBiałaczka.pqs)

2014/08/22 20:00

Weryfikacja modelu

Istotność statystyczna poszczególnych zmiennych w modelu (istotność ilorazu hazardów)

Na podstawie współczynnika oraz jego błędu szacunku możemy wnioskować czy zmienna niezależna, dla której ten współczynnik został oszacowany wywiera istotny wpływ na zmienną zależną. W tym celu posługujemy się testem Walda.

Hipotezy:

$\begin{array}{cc} \mathcal{H}_0: & \beta_i=0,\\ \mathcal{H}_1: & \beta_i\ne 0. \end{array}$

lub równoważnie:

$\begin{array}{cc} \mathcal{H}_0: & HR_i=1,\\ \mathcal{H}_1: & HR_i\ne 1. \end{array}$

Statystykę testową testu Walda wyliczamy według wzoru:

$\begin{displaymath} \chi^2=\left(\frac{b_i}{SE_{b_i}}\right)^2 \end{displaymath}$

Statystyka ta ma asymptotycznie (dla dużych liczności) rozkład chi-kwadrat z $1$ stopniem swobody .

Wyznaczoną na podstawie statystyki testowej wartość $p$ porównujemy z poziomem istotności $\alpha$ :

$\begin{array}{ccl} $ jeżeli $ p \le \alpha & \Longrightarrow & $ odrzucamy $ \mathcal{H}_0 $ przyjmując $ \mathcal{H}_1, \\ $ jeżeli $ p > \alpha & \Longrightarrow & $ nie ma podstaw, aby odrzucić $ \mathcal{H}_0. \\ \end{array}$

Jakość zbudowanego modelu

Dobry model powinien spełniać dwa podstawowe warunki: powinien być dobrze dopasowany i możliwie jak najprostszy. Jakość modelu proporcjonalnego hazardu Cox'a możemy ocenić kilkoma ogólnymi miarami, które opierają się na:

$L_{FM}$ - maksimum funkcji wiarygodności modelu pełnego (z wszystkimi zmiennymi),

$L_0$ - maksimum funkcji wiarygodności modelu zawierającego jedynie wyraz wolny,

$d$ - obserwowanej liczbie niepożądanych zdarzeń.

Kryteria informacyjne opierają się na entropii informacji niesionej przez model (niepewności modelu) tzn. szacują utraconą informację, gdy dany model jest używany do opisu badanego zjawiska. Powinniśmy zatem wybierać model o minimalnej wartości danego kryterium informacyjnego.

$AIC$ , $AICc$ i $BIC$ jest rodzajem kompromisu pomiędzy dobrocią dopasowania i złożonością. Drugi element sumy we wzorach na kryteria informacyjne (tzw. funkcja straty lub kary) mierzy prostotę modelu. Zależy on od liczby parametrów w modelu ( $k$ ) i liczby obserwacji kompletnych ( $d$ ). W obu przypadkach element ten rośnie wraz ze wzrostem liczby parametrów i wzrost ten jest tym szybszy im mniejsza jest liczba obserwacji.

Kryterium informacyjne nie jest jednak miarą absolutną, tzn. jeśli wszystkie porównywane modele źle opisują rzeczywistość w kryterium informacyjnym nie ma sensu szukać ostrzeżenia.

Kryterium informacyjne Akaikego (ang. Akaike information criterion)

$\begin{displaymath} AIC=-2\ln L_{FM}+2k, \end{displaymath}$

Jest to kryterium asymptotyczne - odpowiednie dla dużych prób.

Poprawione kryterium informacyjne Akaikego

$\begin{displaymath} AICc=AIC+\frac{2k(k+1)}{d-k-1}, \end{displaymath}$

Poprawka kryterium Akaikego dotyczy wielkości próby (liczby zdarzeń niepożądanych), przez co jest to miara rekomendowana również dla prób o małych licznościach.

Bayesowskie kryterium informacyjne Schwartza (ang. Bayes Information Criterion lub Schwarz criterion)

$\begin{displaymath} BIC=-2\ln L_{FM}+k\ln(d), \end{displaymath}$

Podobnie jak poprawione kryterium Akaikego uwzględnia wielkość próby (liczbę zdarzeń niepożądanych) - Volinsky i Raftery 2000r⁴⁾

Pseudo R $^2$ - tzw. McFadden $R^2$ jest miarą dopasowania modelu (odpowiednikiem współczynnika determinacji wielorakiej $R^2$ wyznaczanego dla liniowej regresji wielorakiej).

Wartość tego współczynnika mieści się w przedziale $<0; 1)$ , gdzie wartości bliskie 1 oznaczają doskonałe dopasowanie modelu, $0$ - zupełny bark dopasowania. Współczynnik $R^2_{Pseudo}$ wyliczamy z wzoru:

$\begin{displaymath} R^2_{Pseudo}=1-\frac{\ln L_{FM}}{\ln L_0}. \end{displaymath}$

Ponieważ współczynnik $R^2_{Pseudo}$ nie przyjmuje wartości 1 i jest wrażliwy na ilość zmiennych w modelu, wyznacza się jego poprawioną wartość:

$\begin{displaymath} R^2_{Nagelkerke}=\frac{1-e^{-(2/d)(\ln L_{FM}-\ln L_0)}}{1-e^{(2/d)\ln L_0}} \quad \textrm{lub}\quad R^2_{Cox-Snell}=1-e^{\frac{(-2\ln L_0)-(-2\ln L_{FM})}{d}}. \end{displaymath}$

Istotność statystyczna wszystkich zmiennych w modelu

Podstawowym narzędziem szacującym istotność wszystkich zmiennych w modelu jest test ilorazu wiarygodności. Test ten weryfikuje hipotezę:

$\begin{array}{cc} \mathcal{H}_0: & \textrm{wszystkie }\beta_i=0,\\ \mathcal{H}_1: & \textrm{istnieje }\beta_i\neq0. \end{array}$

Statystyka testowa ma postać:

$\begin{displaymath} \chi^2=-2\ln(L_0/L_{FM})=-2\ln(L_0)-(-2\ln(L_{FM})). \end{displaymath}$

Statystyka ta ma asymptotycznie (dla dużych liczności) rozkład chi-kwadrat z $k$ stopniami swobody.

Wyznaczoną na podstawie statystyki testowej wartość $p$ porównujemy z poziomem istotności $\alpha$ :

$\begin{array}{ccl} $ jeżeli $ p \le \alpha & \Longrightarrow & $ odrzucamy $ \mathcal{H}_0 $ przyjmując $ \mathcal{H}_1, \\ $ jeżeli $ p > \alpha & \Longrightarrow & $ nie ma podstaw, aby odrzucić $ \mathcal{H}_0. \\ \end{array}$

AUC – pole pod krzywą ROC – Krzywa ROC – zbudowana na podstawie informacji o wystąpieniu zdarzenia lub jego braku oraz kombinacji zmiennych niezależnych i parametrów modelu – pozwala na ocenę zdolności zbudowanego modelu regresji Cox'a do klasyfikacji przypadków do dwóch grup: (1 – zdarzenie) i (0 – brak zdarzenia). Powstała w ten sposób krzywa, a w szczególności pole pod nią, obrazuje jakość klasyfikacyjną modelu. Gdy krzywa ROC pokrywa się z przekątną $y = x$ , to decyzja o przyporządkowaniu przypadku do wybranej klasy (1) lub (0) podejmowana na podstawie modelu jest tak samo dobra jak losowy podział badanych przypadków do tych grup. Jakość klasyfikacyjna modelu jest dobra, gdy krzywa znajduje się znacznie powyżej przekątnej $y=x$ , czyli gdy pole pod krzywą ROC jest znacznie większe niż pole pod prostą $y=x$ , zatem większe niż $0.5$

Hipotezy:

$\begin{array}{cl} \mathcal{H}_0: & AUC=0.5, \\ \mathcal{H}_1: & AUC\neq 0.5. \end{array}$

Statystyka testowa ma postać:

$\begin{displaymath} Z=\frac{AUC-0.5}{SE_{0.5}}, \end{displaymath}$ gdzie:

$SE_{0.5}$ - błąd pola.

Statystyka $Z$ ma asymptotycznie (dla dużych liczności) rozkład normalny.

Wyznaczoną na podstawie statystyki testowej wartość $p$ porównujemy z poziomem istotności $\alpha$ :

$\begin{array}{ccl} $ jeżeli $ p \le \alpha & \Longrightarrow & $ odrzucamy $ \mathcal{H}_0 $ przyjmując $ \mathcal{H}_1, \\ $ jeżeli $ p > \alpha & \Longrightarrow & $ nie ma podstaw, aby odrzucić $ \mathcal{H}_0. \\ \end{array}$

Dodatkowo, dla krzywej ROC podawana jest proponowana wartość punktu odcięcia kombinacji zmiennych niezależnych i parametrów modelu.

Przykład c.d. (plik: remisjaBiałaczka.pqs)

2014/08/22 20:00

Analiza reszt modelu

Analiza reszt modelu pozwala na weryfikację jego założeń. Głównym jej celem w regresji Cox'a jest lokalizacja wartości odstających i badanie proporcjonalności hazardu. Standardowo w modelach regresji reszty oblicza się jako różnice obserwowanych i przewidywanych przez model wartości zmiennej zależnej. Jednakże w przypadku występowania zmiennych uciętych taka koncepcja wyznaczania reszt nie jest odpowiednia. W programie można analizować reszty opisywane jako: Martingale, Deviance, Schoenfeld. Reszty te można wyrysowywać względem czasu lub zmiennych niezależnych.

Założenie proporcjonalności hazardu

Opracowano szereg graficznych metod pozwalających na ocenę adekwatności modelu proporcjonalnego hazardu (Lee i Wang 2003⁵⁾). Najszerzej stosowane są metody bazujące na resztach modelu. Podobnie jak inne graficzne metody oceny proporcjonalności hazardu jest to metoda subiektywna. Aby założenie proporcjonalnego hazardu było spełnione reszty względem czasu nie powinny układać się w żaden wzór ale powinny być losowo rozłożone w okolicach wartości 0.

Martingale - reszty te mogą być interpretowane jako różnica w czasie $[0,t]$ pomiędzy obserwowaną licznością zdarzeń niepożądanych i przewidywaną przez model licznością tych zdarzeń. Wartość oczekiwana dla tych reszt to 0, ale ich wadą jest skośny rozkład utrudniający interpretację ich wykresu (zawierają się w przedziale od $-\infty$ do 1).

Deviance - podobnie jak Martingale asymptotycznie uzyskują wartość oczekiwaną równą 0, ale rozkładają się symetrycznie wokół zera z odchyleniem standardowym równym 1, gdy dopasowywany model jest odpowiedni. Wartość Deviance jest dodatnia, gdy obiekt badany przeżywa krócej niż oczekiwany na podstawie modelu czas, a ujemna gdy ten czas jest dłuższy. Analiza tych reszt jest wykorzystywana w badaniu proporcjonalności hazardu - ale jest to przede wszystkim narzędzie identyfikujące wartości odstające. W raporcie reszt, te z nich, które są oddalone o więcej niż 3 odchylenia standardowe od 0 oznaczone są kolorem czerwonym.

Schoenfeld - reszty te wyliczane są oddzielnie dla każdej zmiennej niezależnej i definiowane tylko dla obserwacji kompletnych. Dla każdej zmiennej niezależnej suma reszt Schoenfeld'a i ich wartość oczekiwana to 0. Zaletą przedstawiania reszt względem czasu dla każdej zmiennej jest możliwość zidentyfikowania zmiennej nie spełniającej w modelu założenia proporcjonalności hazardu. Jest to ta zmienna, której wykres reszt układa się w systematyczny wzór (najczęściej bada się tu liniową zależność reszt od czasu). Równomierne rozłożenie punktów względem wartości 0 świadczy o braku zależności reszt od czasu - czyli spełnieniu założenia proporcjonalności hazardu przez daną zmienną w modelu.

Gdy dla którejkolwiek zmiennej w modelu Cox'a nie jest spełnione założenie proporcjonalności hazardu, jednym z możliwych rozwiązań jest wykonanie analizy Cox'a oddzielnie dla każdego poziomu tej zmiennej.

Przykład c.d. (plik: remisjaBiałaczka.pqs)

2014/08/22 20:00

¹⁾

Cox D.R. (1972), Regression models and life tables. Journal of the Royal Statistical Society, B34:187-220

²⁾

Breslow N.E. (1974), Covariance analysis of censored survival data. Biometrics, 30(1):89–99

³⁾

Peduzzi P., Concato J., Feinstein A.R., Holford T.R. (1995), Importance of events per independent variable in proportional hazards regression analysis. II. Accuracy and precision of regression estimates. Journal of Clinical Epidemiology, 48:1503-1510

⁴⁾

Volinsky C.T., Raftery A.E. (2000) , Bayesian information criterion for censored survival models. Biometrics, 56(1):256–262

⁵⁾

Lee E. T., Wang J. W. (2003), Statistical Methods for Survival Data Analysis, ed. third, Wiley

PQStat - Baza Wiedzy

Narzędzia użytkownika

Narzędzia witryny

Pasek boczny

Spis treści

Regresja proporcjonalnego hazardu Cox'a

Iloraz Hazardu (HR)

Weryfikacja modelu

Analiza reszt modelu

Narzędzia strony