Regresja logistyczna

Kiedy stosować:

Budowany model regresji logistycznej (podobnie jak liniowej regresji wielorakiej) pozwala na zbadanie wpływu wielu zmiennych niezależnych (\(X_1, X_2,..., X_k\)) na jedną zmienną zależną (\(Y\)). Tym razem jednak zmienna zależna przyjmuje jedynie dwie wartości, np. chory/zdrowy, niewypłacalny/wypłacalny itp. Owe dwie wartości kodowane są jako (1)/(0), gdzie (1) to wartość wyróżniona – posiadanie danej cechy, a (0) – brak danej cechy.

Funkcja, na której oparty jest model regresji logistycznej, wylicza nie dwupoziomową zmienną \(Y\), a prawdopodobieństwo przyjęcia przez tę zmienną wyróżnionej wartości:

\[ P(Y=1|X_1,X_2,...,X_k)=\frac{e^Z}{1+e^Z}, \]

gdzie:
\(P(Y=1|X_1,X_2,...,X_k)\) – prawdopodobieństwo przyjęcia wartości wyróżnionej (1) pod warunkiem uzyskania konkretnych wartości zmiennych niezależnych, tzw. prawdopodobieństwo przewidywane dla 1,
\(Z\) najczęściej wyrażone jest zależnością liniową: \(Z=\beta_0+\sum_{i=1}^k\beta_iX_i\),
\(X_1,X_2,\ldots X_k\) – zmienne niezależne, objaśniające,
\(\beta_0,\beta_1,\beta_2,\ldots \beta_k\) – parametry.

Logitem nazywamy przekształcenie tego modelu do postaci:

\[ \ln\left(\frac{P}{1-P}\right)=Z. \]

Warunki stosowania

zmienna zależna \(Y\) dychotomiczna, zakodowana jako (1)/(0),
zmienne niezależne niewspółliniowe,
odpowiednio duża liczba obserwacji.

Uwaga! Funkcja \(Z\) może być również opisana zależnością wyższego stopnia, np. kwadratową – do modelu wprowadzamy wówczas zmienną zawierającą kwadrat danej zmiennej niezależnej \(X_i^2\).
Budując model, należy pamiętać, że liczba obserwacji powinna być dziesięciokrotnie większa lub równa liczbie szacowanych parametrów modelu (\(n\ge 10(k+1)\)). Jednak coraz częściej stosuje się bardziej restrykcyjne kryterium zaproponowane przez P. Peduzzi i innych w roku 1996[133], mówiące, iż liczba obserwacji powinna być dziesięciokrotnie większa lub równa stosunkowi liczby zmiennych niezależnych (\(v\)) i mniejszej z proporcji liczności (\(p\)) opisanych przez zmienną zależną (tzn. proporcji chorych lub zdrowych), czyli \(n\ge 10v/p\).
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. Zmienne współliniowe należy usunąć z modelu bądź zbudować z nich jedną zmienną niezależną, np. zamiast współliniowych zmiennych: wiek matki i wiek ojca, można zbudować zmienną wiek rodziców.
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ą liczbę 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).

Definicje miar

Zmienne fikcyjne i interakcje w modelu – 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 regresji – rozwiązaniem równania modelu jest wektor ocen parametrów \(\beta_0,\beta_1,\ldots,\beta_k\), nazywanych współczynnikami regresji (\(b_0, b_1,...,b_k\)), szacowanych poprzez metodę największej wiarygodności, czyli poprzez poszukiwanie maksimum funkcji wiarygodności \(L\) (w programie użyto algorytmu iteracyjnego Newton-Raphson). Na podstawie tych wartości możemy wnioskować o wielkości wpływu zmiennej niezależnej (dla której ten współczynnik został oszacowany) na zmienną zależną. Każdy współczynnik obarczony jest pewnym błędem szacunku, wyliczanym ze wzoru \(SE_b=\sqrt{diag(H^{-1})_b}\), gdzie \(diag(H^{-1})\) to główna przekątna macierzy kowariancji.

Jednostkowy Iloraz Szans – na podstawie współczynników, dla każdej zmiennej niezależnej w modelu, wylicza się łatwą w interpretacji miarę, jaką jest jednostkowy Iloraz Szans:

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

Otrzymany Iloraz Szans wyraża zmianę szansy na wystąpienie wyróżnionej wartości (1), 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ść OR interpretujemy następująco:

\(OR >1\) oznacza stymulujący wpływ badanej zmiennej niezależnej na uzyskanie wyróżnionej wartości (1), tj. mówi, o ile wzrasta szansa na wystąpienie wyróżnionej wartości (1), gdy zmienna niezależna wzrasta o jeden poziom;
\(OR <1\) oznacza destymulujący wpływ badanej zmiennej niezależnej na uzyskanie wyróżnionej wartości (1), tj. mówi, o ile spada szansa na wystąpienie wyróżnionej wartości (1), gdy zmienna niezależna wzrasta o jeden poziom;
\(OR\approx1\) oznacza, że badana zmienna niezależna nie ma wpływu na uzyskanie wyróżnionej wartości (1).

U-smile – wykres czterech podklas (nie-zdarzenia: lepsza i gorsza predykcja; zdarzenia: gorsza i lepsza predykcja) z wybranym współczynnikiem (\(rLR\), \(BA\) lub \(RB\)) na osi Y; rozmiar punktu odpowiada udziałowi podklasy \(I\), linie klas są ciągłe przy istotnej poprawie (dla \(rLR\)), a pozioma linia odniesienia oznacza brak zmiany. Uśmiech oznacza poprawę predykcji względem modelu zerowego w obu klasach (patrz U-smile).

Iloraz Szans – wzór ogólny – program PQStat wylicza jednostkowy Iloraz Szans. Jego modyfikacja, na podstawie ogólnego wzoru, umożliwia zmianę interpretacji uzyskanego wyniku. Iloraz szans na wystąpienie stanu wyróżnionego w ogólnym przypadku jest wyliczany jako iloraz dwóch szans wyznaczonych dla dwóch kategorii porównywanej zmiennej niezależnej.

Uwaga! Jeśli analizę przeprowadzamy dla modelu innego niż liniowy lub uwzględniamy interakcję, wówczas na podstawie ogólnego wzoru możemy wyliczyć odpowiedni Iloraz Szans, zmieniając formułę wyrażającą \(Z\).

Jak ustawić analizę

Statystyka zaawansowana \(\to\) Modele wielowymiarowe \(\to\) Regresja logistyczna.

Weryfikacja modelu

Istotność statystyczna poszczególnych zmiennych w modelu (istotność ilorazu szans) – 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.

Jakość zbudowanego modelu – dobry model powinien spełniać dwa podstawowe warunki: powinien być dobrze dopasowany i możliwie jak najprostszy. Jakość modelu regresji logistycznej możemy ocenić kilkoma miarami, opartymi 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, \(n\) – liczności próby:

Kryteria informacyjne (\(AIC\), \(AICc\), \(BIC\)) – opierają się na entropii informacji niesionej przez model. Powinniśmy wybierać model o minimalnej wartości danego kryterium informacyjnego.
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 \(\langle0; 1)\), gdzie wartości bliskie 1 oznaczają doskonałe dopasowanie modelu. Ponieważ współczynnik \(R^2_{Pseudo}\) nie przyjmuje wartości 1 i jest wrażliwy na liczbę zmiennych w modelu, wyznacza się jego poprawione wartości: \(R^2_{Nagelkerke}\) i \(R^2_{Cox-Snell}\).
Istotność statystyczna wszystkich zmiennych w modelu – podstawowym narzędziem szacującym istotność wszystkich zmiennych w modelu jest test ilorazu wiarygodności.
Test Hosmera-Lemeshowa – dla różnych podgrup danych porównuje obserwowane liczności występowania wartości wyróżnionej \(O_g\) i przewidywane prawdopodobieństwo \(E_g\). Jeśli \(O_g\) i \(E_g\) są wystarczająco bliskie, można założyć, że zbudowano dobrze dopasowany model. Do obliczeń najpierw obserwacje są dzielone na \(G\) podgrup, zwykle na decyle (\(G=10\)).
AUC – pole pod krzywą ROC – krzywa ROC, zbudowana w oparciu o wartość zmiennej zależnej oraz przewidywane prawdopodobieństwo zmiennej zależnej \(P\), pozwala na ocenę zdolności zbudowanego modelu regresji logistycznej do klasyfikacji przypadków do dwóch grup: (1) i (0). Dodatkowo, dla krzywej ROC podawana jest proponowana wartość punktu odcięcia prawdopodobieństwa przewidywanego oraz tabela podająca wielkość czułości i swoistości dla każdego możliwego punktu odcięcia.

Uwaga! Więcej możliwości w wyliczeniu punktu odcięcia daje moduł Krzywa ROC. Analizę przeprowadzamy na podstawie wartości obserwowanych i prawdopodobieństwa przewidywanego, które uzyskujemy w analizie regresji logistycznej.

Klasyfikacja – na podstawie wybranego punktu odcięcia prawdopodobieństwa przewidywanego można sprawdzić jakość klasyfikacji. Punkt odcięcia to domyślnie wartość 0.5. Użytkownik może zmienić tę wartość na dowolną wartość z przedziału \((0,1)\), np. wartość sugerowaną przez krzywą ROC. W wyniku uzyskamy tabelę klasyfikacji oraz procent poprawnie zaklasyfikowanych przypadków, procent poprawnie zaklasyfikowanych (0) – swoistość, oraz procent poprawnie zaklasyfikowanych (1) – czułość.

U-smile – ocena modelu względem modelu zerowego – metoda U-smile (Universal Smile Layout for Explanation, Kubiak i in. 2024[178], Więckowska i in. 2025[179], 2026[180]) porównuje dwa modele przewidujące prawdopodobieństwo zdarzenia na tych samych obserwacjach. W oknie regresji logistycznej modelem odniesienia jest model zerowy (tylko wyraz wolny, czyli częstość zdarzeń w próbie), a modelem nowym – zbudowany model. Obserwacje dzieli się na cztery podklasy:

nie-zdarzenia (0) z lepszą predykcją (prawdopodobieństwo zmalało),
nie-zdarzenia (0) z gorszą predykcją,
zdarzenia (1) z gorszą predykcją,
zdarzenia (1) z lepszą predykcją (prawdopodobieństwo wzrosło).

Dla każdej podklasy wyznaczany jest:

współczynnik jakościowy dotyczący liczby zmian \(I\) (udział obserwacji podklasy w klasie),
oraz do wyboru jeden ze współczynników ilościowych dotyczących wielkości zmian:
- \(rLR\) (względny iloraz wiarygodności) – znormalizowany do przedziału \(\langle 0;1\rangle\),
- \(BA\) (Brier Alteration) – średnia zmiana kwadratu reszty,
- \(RB\) (Relative Brier) – zmiana względem błędu Briera modelu odniesienia.

Współczynniki netto dla klas (poprawa minus pogorszenie) i współczynnik ogólny pokazują, o ile model poprawia predykcję osobno dla zdarzeń i nie-zdarzeń.
Istotność sprawdza test ilorazu wiarygodności (ogółem), a dla \(rLR\) dodatkowo testy dla każdej klasy oparte na rozkładzie gamma.
Wzory: dla zainteresowanych.

Wykresy w regresji logistycznej

Iloraz szans \(\pm\) przedział ufności – wykres przedstawiający OR wraz z 95% przedziałem ufności dla wyniku każdej zmiennej zwróconej w zbudowanym modelu. Dla zmiennych kategorialnych linia na poziomie 1 wskazuje wartość ilorazu szans dla kategorii referencyjnej.
Wartości obserwowane / Prawdopodobieństwo oczekiwane – wykres przedstawiający wyniki przewidywanego dla każdej osoby prawdopodobieństwa wystąpienia zdarzenia (oś X) oraz wartości prawdziwej, czyli wystąpienia zdarzenia (wartość 1 na osi Y) lub braku zdarzenia (wartość 0 na osi Y). Jeśli model bardzo dobrze prognozuje, to przy lewej stronie wykresu punkty będą się kumulowały w dolnej części, a przy prawej stronie – w górnej.
Krzywa ROC – wykres zbudowany w oparciu o wartość zmiennej zależnej oraz przewidywane prawdopodobieństwo wystąpienia zdarzenia.
Wykres reszt Pearsona – wykres pozwalający ocenić, czy występują odstające dane. Reszty to różnice wartości obserwowanej i przewidywanego przez model prawdopodobieństwa. Wykresy surowych reszt z regresji logistycznej trudno interpretować, dlatego ujednolica się je, wyznaczając reszty Pearsona (surowa reszta podzielona przez pierwiastek kwadratowy z funkcji wariancji). Znak (dodatni lub ujemny) wskazuje, czy obserwowana wartość jest wyższa czy niższa niż wartość dopasowana do modelu, a wielkość wskazuje stopień odchylenia. Reszty Pearsona mniejsze niż -3 lub większe niż 3 sugerują zbyt duże odchylenie danego obiektu.
Jednostkowe zmiany ilorazu szans – wykres przedstawiający serie ilorazów szans wraz z przedziałem ufności, wyznaczane dla każdego z możliwych punktów odcięcia zmiennej umieszczonej na osi X. Umożliwia użytkownikowi wybór jednego, dobrego punktu odcięcia i następnie zbudowanie na tej podstawie nowej dwuwartościowej zmiennej. Wykres dedykowany jest dla oceny zmiennych ciągłych w analizie jednoczynnikowej.
Profile ilorazu szans – wykres przedstawiający serie ilorazów szans wraz z przedziałem ufności, wyznaczane dla wskazanej wielkości okna, tzn. porównujące częstości wewnątrz okna z częstościami umieszczonymi na zewnątrz okna. Najlepiej sprawdza się, gdy poszukujemy funkcji U-kształtnej, tzn. o wysokim ryzyku przy niskich i przy wysokich wartościach badanej zmiennej, a o niskim – przy wartościach przeciętnych. Im szersze okno, tym większe uogólnienie wyników i gładsza funkcja ilorazu szans; im węższe okno, tym bardziej szczegółowe, ale bardziej rozchwiane wyniki. Wykres dedykowany jest dla oceny zmiennych ciągłych w analizie jednoczynnikowej.
Wykres U-smile – na wykresie U-smile cztery podklasy tworzą „uśmiech”, gdy model poprawia predykcję w obu klasach:
- (oś Y) wielkość uśmiechu pokazuje rozmiar poprawy (do wyboru współczynniki \(rLR\), \(BA\), \(RB\)),
- (oś X) od lewej strony przedstawia klasę braku zdarzeń (kolor niebieski), od prawej klasę zdarzeń (kolor czerwony); punkty skrajne to podklasa: poprawa predykcji, a punkty środkowe to podklasa: pogorszenie predykcji,
- symetria (lewa strona wykresu podobna do prawej) oznacza, że poprawa/pogorszenie predykcji dotyczy obu klas w taki sam sposób,
- wielkość punktu odpowiada współczynnikowi \(I\), czyli wskazuje, jak duży odsetek obiektów znajduje się w każdej podklasie,
- linia klasy jest ciągła, gdy poprawa w tej klasie jest istotna (tylko \(rLR\)).
Przykłady możliwych kształtów wykresu U-smile. Poprawa przewidywania: (a) dla obu klas wyników, (b) tylko dla braku zdarzeń, (c) tylko dla zdarzeń. Pogorszenie przewidywania: (d) dla obu klas wyników, (e) tylko dla braku zdarzeń, (f) tylko dla zdarzeń. Okręgi leżące na mniej więcej stałym poziomie oznaczają brak poprawy lub pogorszenia przewidywania w porównaniu z modelem odniesienia (g). Zygzak wskazuje na poprawę przewidywania w jednej klasie wyników i pogorszenie przewidywania w drugiej (h i i).

Wykonanie i interpretacja

Przykład 1. (plik Profile OR.pqs)
Badamy ryzyko występowania choroby A i choroby B w zależności od BMI pacjenta. Ponieważ BMI jest zmienną ciągłą, jej umieszczenie w modelu skutkuje wyznaczeniem jednostkowego ilorazu szans wyznaczającego liniowy trend wzrostu lub spadku ryzyka. Nie wiemy, czy model liniowy będzie dobrym modelem dla analizy tego ryzyka, dlatego przed budowaniem wielowymiarowych modeli regresji logistycznej zbudujemy kilka modeli jednowymiarowych prezentujących tę zmienną na wykresach, by móc ocenić kształt badanej zależności. Do tego celu posłużą wykresy jednostkowych zmian ilorazu szans i profili ilorazu szans, przy czym dla profili wybierzemy okno o wielkości 100, ponieważ prawie każdy pacjent ma inne BMI, więc około 100 pacjentów znajdzie się w każdym oknie.

Choroba A
Jednostkowe zmiany ilorazu szans pokazują, że gdy punkt odcięcia BMI wybierzemy gdzieś między 27 a 37, to uzyskamy istotny statystycznie i dodatni iloraz szans, pokazujący że osoby mające BMI powyżej tej wartości mają istotnie wyższe ryzyko choroby niż osoby poniżej tej wartości.

Profile ilorazu szans pokazują, że czerwona krzywa znajduje się wciąż blisko jedynki, nieco wyżej jest tylko końcówka krzywej, co wskazuje, że może być trudno podzielić BMI na więcej niż 2 kategorie i wybrać dobrą kategorię referencyjną.

podział na 2 wartości: \(OR[95\%CI]=[1.41, 4.90]\), \(p=0.0024\); jednostkowy OR: \(OR[95\%CI]=1.07[1.02, 1.13]\), \(p=0.0052\)

Osoby z BMI powyżej 30 mają istotnie wyższe ryzyko choroby A; alternatywnie można pozostać przy jednostkowym ilorazie szans, wskazującym stały wzrost ryzyka przy wzroście BMI o jednostkę.

Choroba B
Jednostkowe zmiany ilorazu szans pokazują, że gdy punkt odcięcia BMI wybierzemy gdzieś między 22 a 35, to uzyskamy istotny statystycznie i dodatni iloraz szans, pokazujący że osoby mające BMI powyżej tej wartości mają istotnie wyższe ryzyko choroby niż osoby poniżej tej wartości.

Profile ilorazu szans pokazują, że znacznie lepiej byłoby podzielić BMI na 2 lub 4 kategorie, przy czym kategorią referencyjną powinna być kategoria obejmująca BMI gdzieś pomiędzy 19 a 25, ponieważ to ta kategoria znajduje się najniżej i jest mocno oddalona od wyników dla BMI znajdujących się na lewo i na prawo od tego przedziału. Widzimy wyraźny kształt przypominający literę U, co oznacza, że ryzyko choroby jest wysokie przy niskim i przy wysokim BMI.

Podsumowując, mimo że zależność dla jednostkowego ilorazu szans (zależność liniowa) jest istotna statystycznie, nie warto budować takiego właśnie modelu. Znacznie lepiej podzielić BMI na kategorie, gdzie wartością odniesienia będzie przeciętne BMI. Wykorzystując standardowy podział BMI i ustanawiając kategorią odniesienia BMI w normie, uzyskamy ponad 15-krotnie wyższe ryzyko dla osób z niedowagą (OR[95%CI]=15.14[6.93, 33.10]), ponad dziesięciokrotnie dla osób z nadwagą (OR[95%CI]=10.35[6.74, 15.90]) i ponad dwunastokrotnie dla osób z otyłością (OR[95%CI]=12.22[6.94, 21.49]).

Na wykresie ilorazów szans norma BMI wskazana jest na poziomie 1, jako kategoria referencyjna. Dorysowaliśmy linie łączące uzyskane OR i również normę, tak by pokazać, że uzyskany kształt zależności jest tożsamy z wyznaczonym wcześniej poprzez profil ilorazu szans.

Przykład 2. (plik wada.pqs)
Przeprowadzono badanie mające na celu identyfikację czynników ryzyka pewnej rzadko występującej wady wrodzonej u dzieci. W badaniu wzięło udział 395 matek dzieci z tą wadą oraz 375 matek dzieci zdrowych. Zebrane dane to: miejsce zamieszkania, płeć dziecka, masa urodzeniowa dziecka, wiek matki, kolejność ciąży, przebyte poronienia samoistne, infekcje oddechowe, palenie tytoniu, wykształcenie matki.

Budujemy model regresji logistycznej, by sprawdzić, które zmienne mogą wywierać istotny wpływ na występowanie wady. Jako zmienną zależną ustawiamy kolumnę GRUPA, wartością wyróżnioną w tej zmiennej jako 1 jest grupa „badana”, czyli matki dzieci z wadą wrodzoną. Kolejnych 9 zmiennych to zmienne niezależne: MiejsceZam (2=miasto/1=wieś), Płeć (1=mężczyzna/0=kobieta), MasaUr (w kilogramach, z dokładnością do 0.5 kg), WiekM (w latach), KolCiąży (dziecko z której ciąży), PoronSamo (1=tak/0=nie), InfOddech (1=tak/0=nie), Palenie (1=tak/0=nie), WyksztM (1=podstawowe lub niżej/2=zawodowe/3=średnie/4=wyższe).

\(R^2_{Pseudo}=0.11\), test ilorazu wiarygodności: \(p<0.0001\), test Hosmera-Lemeshowa: \(p=0.2753\)

Jakość dopasowania modelu nie jest wysoka, ale model jest istotny statystycznie, a wynik testu Hosmera-Lemeshowa (brak istotności jest tu pożądany) wskazuje na dobre dopasowanie liczności obserwowanych do przewidywanego prawdopodobieństwa.

Interpretacja poszczególnych zmiennych w modelu zaczyna się od sprawdzenia ich istotności. W tym przypadku zmienne, które w istotny sposób są związane z występowaniem wady, to: Płeć (\(p=0.0063\)), MasaUr (\(p=0.0188\)), KolCiąży (\(p=0.0035\)), InfOddech (\(p<0.0001\)), Palenie (\(p=0.0003\)).

Badana wada wrodzona jest wadą rzadką, ale szansa na jej wystąpienie zależy od wymienionych zmiennych w sposób opisany poprzez iloraz szans:

zmienna Płeć: \(OR[95\%CI]=1.60[1.14;2.22]\) – szansa wystąpienia wady u chłopca jest 1.6-krotnie większa niż u dziewczynki;
zmienna MasaUr: \(OR[95\%CI]=0.74[0.57;0.95]\) – im wyższa masa urodzeniowa, tym szansa wystąpienia wady u dziecka jest mniejsza;
zmienna KolCiąży: \(OR[95\%CI]=1.34[1.10;1.63]\) – szansa wystąpienia wady u dziecka wzrasta wraz z każdą kolejną ciążą 1.34-krotnie;
zmienna InfOddech: \(OR[95\%CI]=4.46[2.59;7.69]\) – szansa wystąpienia wady u dziecka, gdy matka w czasie ciąży przechodziła infekcję oddechową, jest 4.46-krotnie większa niż gdyby jej nie przechodziła;
zmienna Palenie: \(OR[95\%CI]=4.44[1.98;9.96]\) – matka paląca w czasie ciąży zwiększa 4.44-krotnie szansę na wystąpienie wady u dziecka.

W przypadku zmiennych nieistotnych statystycznie przedział ufności dla Ilorazu Szans zawiera jedynkę, co oznacza, że zmienne te nie zwiększają ani nie zmniejszają szansy na wystąpienie badanej wady. Nie można więc interpretować uzyskanego ilorazu w podobny sposób jak dla zmiennych istotnych statystycznie. Wpływ poszczególnych zmiennych niezależnych na występowanie wady możemy również opisać przy pomocy wykresu dotyczącego ilorazu szans:

Kontynuacja przykładu 2 (plik wada.pqs): zbudujemy raz jeszcze model regresji logistycznej, ale tym razem zmienną wykształcenie rozbijemy na zmienne fikcyjne (kodowanie zero-jedynkowe). Tracimy tym samym informację o uporządkowaniu kategorii wykształcenia, ale zyskujemy możliwość wnikliwszej analizy poszczególnych kategorii. Rozbicia na zmienne fikcyjne dokonujemy, wybierając w oknie analizy Zm. fikcyjne:

Wykształcenie podstawowe wybieramy jako kategorię odniesienia.

W rezultacie zmienne opisujące wykształcenie stają się istotne statystycznie. Dopasowanie modelu nie ulega znacznej zmianie, ale zmienia się sposób interpretacji ilorazu szans dla wykształcenia:

Zmienna	\(OR[95\%CI]\)
Wykształcenie podstawowe	kategoria referencyjna
Wykształcenie zawodowe	\(0.51[0.26;0.99]\)
Wykształcenie średnie	\(0.42[0.22;0.80]\)
Wykształcenie wyższe	\(0.45[0.22;0.92]\)

Szansa na wystąpienie badanej wady w każdej kategorii wykształcenia odnoszona jest zawsze do szansy wystąpienia wady przy wykształceniu podstawowym. Widzimy, że dla bardziej wykształconych matek iloraz szans jest niższy. Dla matki z wykształceniem zawodowym szansa wystąpienia wady u dziecka stanowi 0.51 część szansy dotyczącej matki z wykształceniem podstawowym; ze średnim – 0.42 część; z wyższym – 0.45 część.

Przykład 3. (plik zadanie.pqs)
Przeprowadzono eksperyment mający na celu zbadanie umiejętności koncentracji grupy dorosłych podczas sytuacji niekomfortowych. W eksperymencie wzięło udział 190 osób (130 osób to zbiór uczący, 40 osób to zbiór testowy). Każda badana osoba dostała pewne zadanie, którego rozwiązanie wymagało skupienia uwagi. Podczas eksperymentu niektóre osoby zostały poddane działaniu czynnika zakłócającego, jakim była podwyższona temperatura powietrza do 32 stopni Celsjusza. Osoby biorące udział w eksperymencie zapytano dodatkowo o ich miejsce zamieszkania, płeć, wiek i wykształcenie. Czas na rozwiązanie zadania ograniczono do 45 minut. Dla osób, które skończyły przed czasem, odnotowano rzeczywisty czas poświęcony na rozwiązanie. Całość naszych obliczeń wykonamy tylko dla osób należących do zbioru uczącego.

Zmienna ROZWIĄZANIE (tak/nie) zawiera wynik eksperymentu, czyli informację o tym, czy zadanie zostało rozwiązane poprawnie czy też nie. Pozostałe zmienne, które mogły wpływać na wynik eksperymentu, to: MIEJSCEZAM (1=miasto/0=wieś), PŁEĆ (1=kobieta/0=mężczyzna), WIEK (w latach), WYKSZTAŁCENIE (1=podstawowe, 2=zawodowe, 3=średnie, 4=wyższe), CZAS rozwiązywania (w minutach), ZAKŁÓCENIA (1=tak/0=nie). Na bazie wszystkich zmiennych zbudowano model regresji logistycznej, gdzie jako stan wyróżniony zmiennej ROZWIĄZANIE wybrano „tak”.

\(R^2_{Pseudo}=0.27\), \(R^2_{Nagelkerke}=0.41\), \(R^2_{Cox-Snell}=0.30\); test Hosmera-Lemeshowa: \(p=0.1725\); test ilorazu wiarygodności: \(p<0.0001\)

Model jest istotny statystycznie i wystarczająco dobrze dopasowany.

W modelu zmienne, które w sposób istotny wpływają na wynik, to: WIEK (\(p=0.0014\)), CZAS (\(p=0.0012\)), ZAKŁÓCENIA (\(p=0.0001\)). Przy czym im osoba rozwiązująca jest młodsza, czas rozwiązywania krótszy i brak jest czynnika zakłócającego, tym większe prawdopodobieństwo poprawnego rozwiązania: WIEK \(OR[95\%CI]=0.90[0.85;0.96]\), CZAS \(OR[95\%CI]=0.91[0.87;0.97]\), ZAKŁÓCENIA \(OR[95\%CI]=0.15[0.06;0.37]\).

Jeśli model miałby zostać użyty do prognozowania, należy ocenić, czy będzie w tym lepszy niż klasyfikacja losowa. Narysujemy wykres U-smile na bazie współczynnika określającego jakość dopasowania \(rLR\). Wykres się „uśmiecha”. Wiarygodność prognozy wskazującej, że osoba rozwiąże zadanie (prawa część wykresu), jest 0.23 razy większa niż w modelu losowym i netto (tzn. po odjęciu osób, u których wiarygodność spadła) dotyczy 58% osób. Wykres jest symetryczny i podobne wyniki dotyczą lewej strony wykresu, czyli nierozwiązania zadania: \(rLR_0=0.3\), \(I_0=55\%\).

Jakość klasyfikacji ocenimy na bazie krzywej ROC.

\(AUC=0.83\), \(p<0.0001\) (istotnie większe niż 0.5)

Na podstawie zbudowanego modelu można klasyfikować.

Proponowany punkt odcięcia dla krzywej ROC wynosi \(0.6949\) i jest nieco wyższy niż standardowo używany w regresji poziom \(0.5\). Klasyfikacja wyznaczona na bazie tego punktu odcięcia daje \(79.23\%\) przypadków zaklasyfikowanych poprawnie, z czego poprawnie zaklasyfikowanych wartości „tak” jest \(72.73\%\) (czułość), wartości „nie” jest \(88.68\%\) (swoistość). Klasyfikacja uzyskana na podstawie standardowej wartości daje nieco mniej, bo \(73.85\%\) przypadków zaklasyfikowanych poprawnie, ale uzyskamy dzięki niej więcej poprawnie zaklasyfikowanych wartości „tak” (\(83.12\%\)), choć mniej poprawnie zaklasyfikowanych wartości „nie” (\(60.38\%\)).

Na tym etapie możemy zakończyć analizę klasyfikacji, lub jeśli wynik nie jest wystarczający, bardziej wnikliwą analizę krzywej ROC możemy przeprowadzić w module Krzywa ROC.
Jednak pamiętajmy, że wszelkich ocen dokonywaliśmy jedynie na zbiorze uczącym. Jeśli model ma rzeczywiście służyć do analizy danych, na których się nie uczył, powinniśmy go ocenić na niezależnym zbiorze testowym. Zostało to zrobione w module Walidacja i predykcja dla regresji logistycznej, gdzie potwierdzono wyniki zarówno na wykresie U-smile, jak i na krzywej ROC.

Ponieważ uznaliśmy, że klasyfikacja na podstawie modelu jest zadowalająca, możemy wyliczyć prognozowaną wartość zmiennej zależnej dla dowolnie zadanych warunków. Sprawdźmy, jakie szanse na rozwiązanie zadania ma osoba, dla której: MIEJSCEZAM (1=miasto), PŁEĆ (1=kobieta), WIEK (50 lat), WYKSZTAŁCENIE (1=podstawowe), CZAS rozwiązywania (20 minut), ZAKŁÓCENIA (1=tak). W tym celu na podstawie wartości współczynnika \(b\) wyliczane jest prawdopodobieństwo przewidywane:

\[ \begin{array}{l} P(Y=tak|MIEJSCEZAM,PŁEĆ,WIEK,WYKSZTAŁCENIE,CZAS,ZAKŁÓCENIA)=\\[0.2cm] =\frac{e^{7.23-0.45\textrm{\scriptsize \textit{MIEJSCEZAM}}-0.45\textrm{\scriptsize\textit{PŁEĆ}}-0.1\textrm{\scriptsize\textit{WIEK}}+0.46\textrm{\scriptsize\textit{WYKSZTAŁCENIE}}-0.09\textrm{\scriptsize\textit{CZAS}}-1.92\textrm{\scriptsize\textit{ZAKŁÓCENIA}}}}{1+e^{7.23-0.45\textrm{\scriptsize\textit{MIEJSCEZAM}}-0.45\textrm{\scriptsize\textit{PŁEĆ}}-0.1\textrm{\scriptsize\textit{WIEK}}+0.46\textrm{\scriptsize\textit{WYKSZTAŁCENIE}}-0.09\textrm{\scriptsize\textit{CZAS}}-1.92\textrm{\scriptsize\textit{ZAKŁÓCENIA}}}}=\\[0.2cm] =\frac{e^{7.231-0.453\cdot1-0.455\cdot1-0.101\cdot50+0.456\cdot1-0.089\cdot20-1.924\cdot1}}{1+e^{7.231-0.453\cdot1-0.455\cdot1-0.101\cdot50+0.456\cdot1-0.089\cdot20-1.924\cdot1}} \end{array} \]

W rezultacie tych obliczeń program zwróci wynik:

Uzyskane prawdopodobieństwo rozwiązania zadania wynosi \(0.1215\), więc na podstawie punktu odcięcia \(0.5\) przewidziany wynik to \(0\) – czyli zadanie nierozwiązane poprawnie.

DLA ZAINTERESOWANYCH

Iloraz szans – wzór ogólny
Wzory:
Dla zmiennej niezależnej \(X_1\), przy \(Z\) wyrażonym zależnością liniową, wyliczamy szansę dla kategorii pierwszej:

\[ Szansa(1)=\frac{P(1)}{1-P(1)}=e^Z(1)=e^{\beta_0+\beta_1X_1(1)+\beta_2X_2+...+\beta_kX_k}, \]

oraz szansę dla kategorii drugiej:

\[ Szansa(2)=\frac{P(2)}{1-P(2)}=e^Z(2)=e^{\beta_0+\beta_1X_1(2)+\beta_2X_2+...+\beta_kX_k}. \]

Iloraz Szans dla zmiennej \(X_1\) wyraża się wówczas wzorem:

\[ \begin{array}{lll} OR_1(2)/(1) &=&\frac{Szansa(2)}{Szansa(1)}=\frac{e^{\beta_0+\beta_1X_1(2)+\beta_2X_2+...+\beta_kX_k}}{e^{\beta_0+\beta_1X_1(1)+\beta_2X_2+...+\beta_kX_k}}\\ &=& e^{\beta_0+\beta_1X_1(2)+\beta_2X_2+...+\beta_kX_k-[\beta_0+\beta_1X_1(1)+\beta_2X_2+...+\beta_kX_k]}\\ &=& e^{\beta_1X_1(2)-\beta_1X_1(1)}=e^{\beta_1[X_1(2)-X_1(1)]}=\\ &=& \left(e^{\beta_1}\right)^{[X_1(2)-X_1(1)]}. \end{array} \]

Przykład: jeśli zmienną niezależną jest wiek wyrażony w latach, to różnica pomiędzy sąsiadującymi kategoriami wieku, np. 25 lat i 26 lat, wynosi 1 rok \(\left(X_1(2)-X_1(1)=26-25=1\right)\). Wówczas otrzymamy jednostkowy Iloraz Szans \(OR=\left(e^{\beta_1}\right)^1\), który mówi, o ile zmieni się szansa na wystąpienie wyróżnionej wartości, gdy wiek zmieni się o 1 rok. Iloraz szans wyliczony dla niesąsiadujących kategorii, np. 25 i 30 lat, będzie pięcioletnim Ilorazem Szans \(OR=\left(e^{\beta_1}\right)^5\), ponieważ różnica \(X_1(2)-X_1(1)=30-25=5\).

Test Walda – istotność poszczególnych zmiennych
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: & OR_i=1,\\ \mathcal{H}_1: & OR_i\ne 1. \end{array}\)

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

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

Statystyka ta ma asymptotycznie (dla dużych liczności) rozkład \(\chi^2\) z 1 stopniem 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\).

Kryterium informacyjne Akaikego

\[ AIC=-2\ln L_{FM}+2k, \]

jest kryterium asymptotycznym – odpowiednim dla dużych prób.

Poprawione kryterium informacyjne Akaikego

\[ AICc=AIC+\frac{2k(k+1)}{n-k-1}. \]

Bayesowskie kryterium informacyjne Schwarza

\[ BIC=-2\ln L_{FM}+k\ln(n). \]

Pseudo R\(^2\)

\[ R^2_{Pseudo}=1-\frac{\ln L_{FM}}{\ln L_0}. \]

\[ R^2_{Nagelkerke}=\frac{1-e^{-(2/n)(\ln L_{FM}-\ln L_0)}}{1-e^{(2/n)\ln L_0}} \quad \textrm{lub}\quad R^2_{Cox-Snell}=1-e^{\frac{(-2\ln L_0)-(-2\ln L_{FM})}{n}}. \]

Test ilorazu wiarygodności
Hipotezy:

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

Wzory:
Statystyka testowa ma postać:

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

Statystyka ta 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\).

Test Hosmera-Lemeshowa
Hipotezy:

\[ \begin{array}{cl} \mathcal{H}_0: & O_g=E_g \text{ dla wszystkich kategorii,}\\ \mathcal{H}_1: & O_g \neq E_g \text{ dla przynajmniej jednej kategorii.} \end{array} \]

Wzory:
Statystyka testowa ma postać:

\[ H=\sum_{g=1}^G\frac{(O_g-E_g)^2}{E_g(1-\frac{E_g}{N_g})}, \]

gdzie \(N_g\) – liczba obserwacji w grupie \(g\).
Statystyka ta ma asymptotycznie (dla dużych liczności) rozkład \(\chi^2\) z \(G-2\) stopniami swobody.

Test istotności AUC
Interpretacja:
Gdy krzywa ROC pokrywa się z przekątną \(y = x\), to decyzja o przyporządkowaniu przypadku do wybranej klasy podejmowana na podstawie modelu jest tak samo dobra jak losowy podział badanych przypadków do tych grup. Jakość klasyfikacyjna modelu jest dobra, gdy pole pod krzywą ROC jest znacznie większe niż 0.5.

Hipotezy:

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

Wzory:
Statystyka testowa ma postać:

\[ Z=\frac{AUC-0.5}{SE_{0.5}}, \]

gdzie \(SE_{0.5}\) – błąd pola. Statystyka \(Z\) ma asymptotycznie (dla dużych liczności) rozkład normalny.

Współczynniki U-smile
Oznaczenia: \(y_i\in\{0,1\}\) – wynik obserwacji, \(p_i^{ref}\) i \(p_i\) – prawdopodobieństwa zdarzenia z modelu odniesienia i modelu nowego, \(n_0\), \(n_1\) – liczności nie-zdarzeń i zdarzeń, \(n=n_0+n_1\). Podklasy: \(0^{+}\) (nie-zdarzenia, \(p_i<p_i^{ref}\)), \(0^{-}\) (nie-zdarzenia, \(p_i>p_i^{ref}\)), \(1^{-}\) (zdarzenia, \(p_i<p_i^{ref}\)), \(1^{+}\) (zdarzenia, \(p_i>p_i^{ref}\)); obserwacje z \(p_i=p_i^{ref}\) nie należą do żadnej podklasy.

Współczynniki \(I\): \(I_0^{+}=n_0^{+}/n_0\), \(I_0^{-}=n_0^{-}/n_0\), \(I_1^{-}=n_1^{-}/n_1\), \(I_1^{+}=n_1^{+}/n_1\); netto \(I_0=I_0^{+}-I_0^{-}\), \(I_1=I_1^{+}-I_1^{-}\); ogółem \(I=I_0+I_1\).

Współczynniki \(BA\) i \(RB\): dla reszt kwadratowych \(r_i=(y_i-p_i)^2\) i \(r_i^{ref}=(y_i-p_i^{ref})^2\) oraz różnic sum kwadratów w podklasie \(S\): \(\Delta SS_S=\sum_{i\in S} r_i^{ref}-\sum_{i\in S} r_i\)

\[ BA_0^{+}=\frac{\Delta SS_{0^{+}}}{n_0},\quad BA_0^{-}=\frac{-\Delta SS_{0^{-}}}{n_0},\quad BA_1^{-}=\frac{-\Delta SS_{1^{-}}}{n_1},\quad BA_1^{+}=\frac{\Delta SS_{1^{+}}}{n_1}. \]

Współczynniki \(RB\) liczone są tak samo, ale z mianownikami \(SS_0^{ref}=\sum_{y_i=0} r_i^{ref}\) i \(SS_1^{ref}=\sum_{y_i=1} r_i^{ref}\). Netto: \(BA_0=BA_0^{+}-BA_0^{-}\), \(BA_1=BA_1^{+}-BA_1^{-}\); ogółem \(BA=\frac{n_0}{n}BA_0+\frac{n_1}{n}BA_1\) (analogicznie \(RB\)).

Współczynnik \(rLR\): dla składników log-wiarygodności \(\ell_i=y_i\ln p_i+(1-y_i)\ln(1-p_i)\) i \(\ell_i^{ref}\) (analogicznie dla \(p_i^{ref}\)) oraz \(\Delta_S=\sum_{i\in S}\ell_i^{ref}-\sum_{i\in S}\ell_i\):

\[ LR_0^{+}=-2\Delta_{0^{+}},\quad LR_0^{-}=2\Delta_{0^{-}},\quad LR_1^{-}=2\Delta_{1^{-}},\quad LR_1^{+}=-2\Delta_{1^{+}},\quad LR_0=LR_0^{+}-LR_0^{-},\quad LR_1=LR_1^{+}-LR_1^{-}. \]

Wartości te są normalizowane przez maksymalny osiągalny iloraz wiarygodności klasy: \(\max LR_0=-2\sum_{y_i=0}\ell_i^{ref}\) (gdy \(LR_0^{+}\ge LR_0^{-}\); w przeciwnym razie z \(\ell_i\) modelu nowego), analogicznie \(\max LR_1\). Stąd \(rLR_0^{\pm}=LR_0^{\pm}/\max LR_0\), \(rLR_1^{\pm}=LR_1^{\pm}/\max LR_1\), \(rLR_0=LR_0/\max LR_0\), \(rLR_1=LR_1/\max LR_1\) oraz \(rLR=(LR_0+LR_1)/\max LR\), gdzie \(\max LR=-2\sum_i\ell_i^{ref}\) (suma po obserwacjach należących do podklas).

Testy: ogółem – test ilorazu wiarygodności o statystyce \(2\left|\sum_i\ell_i-\sum_i\ell_i^{ref}\right|\), mającej asymptotycznie rozkład \(\chi^2\) z \(df\) równym różnicy liczby parametrów porównywanych modeli (co najmniej 1). Dla \(rLR\) istotność zmiany w klasie sprawdza się z rozkładu gamma: \(|LR_1|\sim\Gamma\left(\frac{df}{2},\frac{1}{2c_1}\right)\) i \(|LR_0|\sim\Gamma\left(\frac{df}{2},\frac{1}{2c_0}\right)\), gdzie \(c_1=n_1/n\), \(c_0=1-c_1\) (parametry kształtu i tempa). Wyznaczone wartości \(p\) porównujemy z poziomem istotności \(\alpha\): gdy \(p \le \alpha\), uznajemy zmianę predykcji za istotną.