Regresja logistyczna dla dowolnej zmiennej niezależnej

Porównując liniowy model prawdopodobieństwa (LMP) z regresją logistyczną (RL), doszliśmy do wniosku, że oba modele staną się tożsame, gdy spełnione będą dwa warunki: występuje tylko jedna zmienna niezależna X oraz zmienna ta jest zero-jedynkowa. Dzieje się tak dlatego, że wtedy oba modele mają tylko 2 punkty w X: 0 i 1, wobec czego naturalnie powstaje zawsze linia prosta (jako funkcja prawdopodobieństwa, P) która je łączy. Poniższy rysunek to pokazuje:

Rysunek nr 1

Już sytuacja, gdy wystąpi wiele zmiennych niezależnych typu 0-1, zmienia ten obraz. Jednak ciekawszy problem powstaje, gdy mowa o dowolnej zmiennej X, nie tylko 0-1. Bo przecież wyprowadzenie RL, które pokazałem tutaj, opierało się w całości na zmiennych binarnych.

Po co RL uogólniać dla dowolnej zmiennej rzeczywistej X? W ekonomii rzadko spotyka się oczywisty podział na tak / nie. Np. moglibyśmy sprawdzić jaka jest szansa, że firma zbankrutuje, jeśli ponosi straty przez 3 lata z rzędu. Ale lepszym sposobem byłaby ocena prawdopodobieństwa bankructwa dla firmy, która ponosi straty dowolną liczbę lat z rzędu. Przypuszczalnie im większa liczba lat ze stratami, tym większa szansa na bankructwo. Innym kryterium mogłaby być wielkość strat: im większe straty, tym większa szansa na upadek firmy. Lub na odwrót: im większe zyski, tym mniejsze zagrożenie upadłością. To kryterium wykorzystam później w przykładzie.

Żeby zrozumieć, dlaczego można podstawić pod X dowolną zmienną, należy powrócić do poprzedniego artykułu, w którym analizowałem przejście od funkcji Y = X do postaci Y = a + bX + e, gdzie a, b to stałe współczynniki, e to składnik losowy. Wtedy p = a + bX. Wiemy, że X równa się 0 lub 1 oraz Y równa się 0 lub 1. Powiedzmy jednak, że chcemy, aby X równało się 0 lub 5, a Y ciągle 0 lub 1. Mielibyśmy dwie możliwości osiągnięcia tego.
I) Zamieniamy współczynnik b na 1/5*b. Wtedy Y = a + 1/5*bX + e, oraz  p = a + 1/5*bX
II) Zamieniamy e na e' = -4/5*bX + e. Wtedy Y = a + bX + e', oraz  p = a + bX + E(e'). Zauważmy, że E(e') = E(-4/5*b*X + e) = -4/5*b*X + E(e) = -4/5*b*X, bo X uznajemy za stałą (jest zawsze znana, 0 lub 5), wartość oczekiwana e wynosi 0. Stąd p =  a + bX + E(e') = a + bX - 4/5b*X = a + 1/5*b*X.

Dostaliśmy to samo w obu sposobach.

Teraz powiedzmy, że są dwie zmienne niezależne X1 i X2, pierwotnie typu 0-1. Ale teraz chcemy, by X1 = 0 lub 5, oraz X2 = 0 lub 10. Wtedy, aby zachować wszystko inne, Y = a + 1/5*b1*X1 + 1/10*b2*X2 + e, oraz p = a + 1/5*b1*X1 + 1/10*b2*X2.

Następnie łączymy obie zmienne w jedną nową X: p = a + b(1/5*b1/b*X1 + 1/10*b2/b*X2) =  a + bX, gdzie X = 1/5*b1/b*X1 + 1/10*b2/b*X2. Powiedzmy, że b1 = 0,5, b2 = 0,7, b = 0,6. Wtedy powstają 3 możliwości rozłożenia X:
Jeśli X1 = 5, a X2 = 10, to X = 1/5*0,5/0,6*5 + 1/10*0,7/0,6*10 = 5/6 + 7/6 = 2.
Jeśli X1 = 5 i X2 = 0, to X = 5/6.
Jeśli X1 = 0, a X2 = 10, to X = 7/6.

Jest też czwarta możliwość, że X1 = 0 i X2 = 0. Najlepiej byłoby, gdyby ta sytuacja była przeznaczona dla Y = 0, a pozostałe dla Y = 1. Oznacza to, że dla p(Y = 1) X powinna mieć tylko 3 różne wartości: 2, 5/6 i 7/6. Wówczas odpowiednio obliczymy prawdopodobieństwo, że Y = 1:
p = a + 0,6*2 = a + 1,2;
p = a + 0,6*5/6 = a + 0,5;
p = a + 0,6*7/6 = a + 0,7.

Możemy z tych 3 wartości p określić przedział dla stałej a. Wiemy, że p < 1, więc na pewno a < -0,2. Ponadto p > 0, więc a > -0,5. W sumie więc a będzie się mieścić między -0,5 a -0,2. Np. a = -0,3. Wtedy:
dla X = 2, p = -0,3 + 1,2 = 0,9,
dla X = 5/6, p = -0,3 + 0,5 = 0,2,
dla X = 7/6, p = -0,3 + 0,7 = 0,4.

Zatem dostajemy dodatnią zależność p od X: dla małej wartości (0,83) p = 0,2, dla większej (1,17) p = 0,4, oraz dla dużej (2) p = 0,9. Sumując, widać, że z początkowych dwóch zmiennych binarnych, doszliśmy do dowolnej zmiennej rzeczywistej X. Jeśli zwiększymy ilość zmiennych do X3, X4..., to liczba możliwych wariantów X będzie się zwiększać.

Zależność funkcyjna P od X będzie wyglądać od teraz tak jak na Rysunku nr 2:

Rysunek nr 2

X może znajdować się w przedziale (a, b). Wszystkie wartości powyżej i poniżej są niedozwolone.

Można zapytać, dlaczego od razu nie zapisałem Y = a + bX + e, gdzie Y to zmienna binarna, a X dowolna rzeczywista, reprezentująca jakąś cechę, np. dochód. Rzecz w tym, że wyprowadzenie LMP w poprzednim artykule zacząłem od tożsamości Y = X, gdzie Y była typu 0-1, a więc X też musiała taka być. Nawet przechodząc do słabszej korelacji Y = a + bX + e, musiałem zachować ten typ. Dopiero więcej zmiennych niezależnych pozwoliło mi na przekształcenie binarnego X w dowolnie rzeczywistą zmienną.

Ponieważ LMP jest funkcją liniową, to dla skrajnych wartości X, p będzie oczywiście zachowywać się mało sensownie. Mając naszą funkcję p = -0,3 + 0,6*X, już dla X = 2,2, dostaniemy p = 1,02 > 1. Albo dla X = 0,4, p = -0,06 < 0. Nie ma powodu, by dochód 2 był możliwy, a już 2,2 nie. Dlatego LMP nie nadaje się dla wielowartościowej zmiennej niezależnej lub dla wielu zmiennych zero-jedynkowych (które prowadzą do tego samego).  

Stąd zamiast LMP zastosujemy RL. Dotychczas mogliśmy RL użyć tylko dla binarnego X. Ale dzięki zastosowaniu powyższego toku rozumowania, dostaniemy model już do bardzo szerokich zastosowań.

Tak więc, mamy regresję logistyczną dla 1 zmiennej niezależnej X typu 0-1:

(1)

Funkcja ta przyjmuje tylko dwie wartości. Dla X = 1:

(2)

Dla X = 0:

(3)

We wzorze (2) i (3) beta(0) znika dlatego, że siedzi w niej x2 = 0 (zob. wyprowadzenie regresji logistycznej tutaj).

Fakt, że występują tylko dwa punkty oznacza, że dostajemy z powrotem Rysunek nr 1, a więc tak jak wcześniej powiedziałem dostajemy de facto LMP.

Dla dwóch zmiennych X1 i X2 typu 0-1 ogólna postać RL będzie następująca:

(4)

i dostaniemy już 4 różne wartości P ze względu na wystąpienie czterech różnych zdarzeń.
Dla X1 = 1 i X2 = 0:

Dla X1 = 0 i X2 = 1:

Dla X1 = 1 i X2 = 1:

Dla X1 = 0 i X2 = 0:

 bo musi zajść:

Podobnie jak dla funkcji jednej zmiennej beta(0) = 0, bo siedzi w niej czwarte zdarzenie, gdy jednocześnie X1 = 0 i X2 = 0.

Czyli widać, że wraz z liczbą zmiennych niezależnych 0-1 wzrasta liczba wariantów P.

Stosując identyczne rozumowanie jak wcześniej dla LMP, zastępujemy teraz zmienne niezależne typu 0-1 jedną dowolną zmienną niezależną. W ten sposób powrócimy do wzoru nr (1), lecz X stanie się już zmienną o dowolnych wartościach rzeczywistych. Dzięki temu linia RL przybierze kształt dystrybuanty rozkładu prawdopodobieństwa, tak jak to pokazuje Rysunek nr 3:

Rysunek nr 3

Zmienna X może od teraz przyjmować dowolną wartość, bo funkcja nie spadnie nigdy poniżej 0 ani nie przekroczy 1.

W tym czasie sam Y ciągle przyjmuje 0 lub 1, tak że zależność od X powinna się prezentować następująco:

Rysunek nr 4

Poziom a oraz b nie są jednoznaczne. Dla a możemy przyjąć, że będzie to wartość, od której zaczynają się straty, tak że firma będzie musiała zbankrutować. Natomiast b możemy uznać za taką wartość od której P zbliża się do 1. Wówczas porównujemy Rysunek nr 4 z Rysunkiem nr 3 oraz Rysunkiem nr 2. Należy zaznaczyć, że Rysunek nr 4 wskazuje jedynie na zależność planowaną, przewidywaną, dlatego w rzeczywistości niektóre punkty przecinające 0 będą przecinać 1, a te które przecinają 1 będą przecinać 0. Potem na przykładzie to zobaczymy.

RL można uogólnić na więcej zmiennych niezależnych. W sumie RL będzie mieć postać:

(5)

gdzie X(k) przyjmuje dowolną wartość rzeczywistą.

Na koniec teoretycznej części warto dodać, że pomiędzy LMP a RL ciągle istnieje pewna zależność. Zapiszmy (5) w prostszej formie:

(6)

Wzór (6) przekształćmy na dwa sposoby a i b:
 a)

b)

Po odjęciu równania (b) od (a) i dalszych prostych przekształceniach otrzymamy model logitowy:

(7)

Ten model jest najczęściej wykorzystywany do oszacowania parametrów RL, zresztą prowadzi do niego metoda największej wiarygodności, za pomocą której znajduje się właśnie parametry.

Model (7) dostarcza cennej informacji. Podstawmy P = 0,5. Dostaniemy ln1 = 0, czyli Z = 0, a to oznacza, że model niczego nie wyjaśnia, tzn. beta(0) = 0 oraz każde X(k) = 0 lub beta(k) = 0. Model praktycznie nie istnieje, więc nie istnieje także przewidywane Y. To prowadzi do wniosku, że przewidywane Y = 0.

A więc doszliśmy do wniosku, że gdy P = 0,5, to przewidywane Y = 0. Jakie to ma praktyczne znaczenie? Otóż prowadzi to do wniosku, że podejmując decyzję, np. o inwestowaniu w daną spółkę, powinniśmy kierować się tym, że spółka nie upadnie z prawdopodobieństwem powyżej 0,5.


Przykład. Część A)
Stworzyłem hipotetyczną listę 98 spółek o rocznych dochodach na akcję, które nie zbankrutowały (oznaczone 1) lub zbankrutowały (0). Ogólnie mamy model ze zmienną zależną Y oznaczoną Czy_bez_upadlosci, która przyjmuje 0 (upadek) lub 1 (bez upadku). Zmienna niezależna X oznaczona jako Dochód_firmy przyjmuje wartości rzeczywiste. Na rysunku nr 4 pokazałem, jak przewidujemy zachowanie Y względem X. Ale wiadomo, że nie zawsze ta zależność będzie spełniona. Popatrzmy więc, jak w przykładzie to wygląda naprawdę:

Na podstawie tego wykresu moglibyśmy spróbować wyznaczyć poziomy a i b, korespondując z Rysunkiem nr 4. Poziom Y = 1 zaczyna się dla X bliskiego 0, stąd X = 0 uznamy na początek za wartość a. Jednocześnie oceniamy, że X = 48 to wartość, po której Y już nie przyjmuje w ogóle 0, zatem 48 może być uznane za wartość b. Jednak nie wszystko się zgadza z Rysunkiem nr 4. Wielu wartościom powyżej X = a przyporządkowane są Y = 0 a nie Y = 1, jak należałoby oczekiwać. Sugeruje to, że poziom a powinien zostać przesunięty. Zobaczymy później czy da się to zrobić. 

Część B)
Na podstawie Przykładu nr 1 zbudujemy LMP.

Y = a +bX +e,
gdzie:
Y = Czy_bez_upadlosci
X = Dochód_firmy
e - składnik losowy.

Linia regresji będzie dana wzorem:

P = a +bX,
gdzie:
P - prawdopodobieństwo, że firma nie upadnie

Stosujemy KMNK:

Wartości prognozy zapisałem jako nową zmienną Pr_Czy_bez_upadlosci_lin. Oto graficzne wyniki prognozy wedle LMP:

Zgodnie z oczekiwaniami im większe dochody, tym szansa, że firma nie zbankrutuje, rośnie. Niestety także dla wielu wartości P < 0 oraz P > 1, co dowodzi, że model jest nieprawidłowy. Gdy firma poniosła duże straty, poniżej -30 zł na akcję, P stało się ujemne. Analogicznie, dla dużych dochodów, powyżej 70 zł, P przekracza 1. Należy więc model skorygować, najlepiej stosując RL.

Część C)
Wykorzystujemy te same dane, aby poprawić prognozy za pomocą RL. Wchodzimy w Model logitowy w gretlu. Pamiętajmy, że mamy tu do czynienia z klasycznym, binarnym logitem, dlatego wchodzimy w Ograniczona zmienna zależna -> model logitowy -> binarny. (Pozostałe warianty przeznaczone są dla zmiennej zależnej nie-zerojedynkowej). Zaznaczamy Dochód_firmy jako regresor, Czy_bez_upadlosci jako zależną, zaznaczamy opcję Pokaż wartość p, i OK.

Otrzymujemy oceny parametrów i inne statystyki:

Model jest istotny statystycznie. Macierz błędów omówiłem w poprzednim artykule; porównując liczbę błędów (sumarycznie 30) z prawidłowymi prognozami (68) oceniamy, że model jest dobry. 

Wartości prognozy zapisałem jako nową zmienną Pr_Czy_bez_upadlosci_log. Graficzne wyniki prognozy wg RL przedstawiają się następująco:

Jest to oczywiście realizacja Rysunku nr 3. RL pięknie dostosowuje duże dochody do zerowego prawdopodobieństwa upadłości oraz duże straty do zagrożenia na poziomie prawie 100%.


Część D)
Powiedzmy, że obserwujemy nową firmę, która w ostatnim roku poniosła stratę 15 zł na akcję. Chcemy za pomocą RL ocenić szanse, że firma przetrwa. W gretlu Dane -> dodaj obserwację -> 1. Potem zaznaczamy Dochód_firmy, potem klikamy Dane -> Edycja wartości -> dopisujemy na koniec -15. Następnie wchodzimy w Prognozę mając ciągle otwarty model logitowy. Patrzymy na ostatnią wartość prognozy:

Czyli jest zaledwie 12% szans, że nowa firma nie upadnie, tzn. na 88% upadnie.

Co by się stało, gdyby ta firma zamiast straty 15 zł, miała 15 zł zysku? Zmieńmy dochód na +15. Teraz prawdopodobieństwo bankructwa wynosi 56%.

Czyli pomimo zysku, firma na ponad 50% i tak upadnie.


Część E)
Wszystko już wiemy, potrafimy modelować prawdopodobieństwo upadłości czy przetrwania, ale musimy teraz odpowiedzieć na pytanie czy inwestować w którąś firmę czy nie. Pierwszym kryterium, które moglibyśmy nazwać koniecznym, będzie wybór takich firm, dla których prawdopodobieństwo przetrwania będzie większe od 0,5. Z jednej strony jest to wybór intuicyjny, z drugiej wynika z samego modelu logitowego, o czym już wcześniej się rozpisałem.

W gretlu do poprzednich danych dopisujemy nową zmienną, np. Czy_bez_upadlosci_hat o formule Pr_Czy_bez_upadlosci_log > 0.5. Ta prosta formuła działa na zasadzie prawda-fałsz, a więc 1-0. Jeżeli zmienna Pr_Czy_bez_upadlosci_log > 0.5, to Czy_bez_upadlosci_hat = 1. W przeciwnym wypadku Czy_bez_upadlosci_hat = 0. To oznacza, że dostaniemy przewidywane Y z Rysunku nr 4 (dopisek hat oznacza z ang. daszek, czyli dostaniemy Y z daszkiem, tj. prognozowane Y właśnie). Tak więc dla zmiennej niezależnej X, zmienna zależna Czy_bez_upadlosci_hat, będzie się zachowywać tak:

  
Czyli dostajemy ten sam obraz co na Rysunku nr 4. Teraz widzimy, że poziom X = a nie wynosi 0, jak początkowo uznaliśmy, ale prawie 20. Czyli minimalnym kryterium, aby w firmę zainwestować, jest osiąganie dochodu 20. Porównując z poprzednim wykresem widać, że dochód = 20 koresponduje z P = 0,5, czyli rzeczywiście te poziomy powyżej 20 zł na akcję posiadają szansę na przetrwanie powyżej 50%.

Oczywiście 50% jest minimalnym poziomem, a więc możemy podwyższyć ten próg. Tym progiem może być taka wartość, dla której X = b, bo od tego poziomu żadna firma nie zbankrutowała. Ten poziom się nie zmieni i ciągle wynosi 48, tak jak wcześniej go wyznaczyliśmy w Części A. Teraz patrzymy, jakie p towarzyszy X = 48. To kryterium moglibyśmy nazwać wystarczającym. Tworzymy nową zmienną Czy_bez_upadlosci_min = Dochod_firmy=48. Formuła ta ponownie działa na zasadzie prawda-fałsz, w ten sposób, że przybierze 1 dla firmy osiągającej dochód 48, a dla pozostałych 0. Najszybciej znajdziemy P zaznaczając  Czy_bez_upadlosci_min oraz Pr_Czy_bez_upadlosci_log, a następnie wyświetlamy wykres rozrzutu XY. Dostałem taki wykres:

W sumie wystarczy, aby firma na ok. 85% była wypłacalna, aby w nią zainwestować.

Liniowy model prawdopodobieństwa vs. Regresja logistyczna

Powiedzmy, że zastanawiamy się czy zainwestować w jakieś akcje. Oczywiście standardowym narzędziem do tego jest wycena zdyskontowanych przepływów pieniężnych lub dywidend. Jeśli model wycenił spółkę na poziomie niższym niż bieżąca kapitalizacja, to dostajemy sygnał kupna. W przeciwnym wypadku dostajemy sygnał sprzedaży. Niestety wycena to problem skomplikowany i oparty na wielu założeniach. Nawet wspomagana przez analizę fundamentalną, w której bierzemy pod uwagę nie tylko czynniki wewnętrzne, ale i otoczenie firmy, może zawieść. Łatwiej zastosować samą AF opartą na wskaźnikach lub analizie dyskryminacyjnej. Zagadnienie tej ostatniej już poruszyłem (tutaj , tutaj , dla banków - tutaj) i należy stwierdzić, że jest wyjątkowo użyteczna, będąc jednocześnie techniką prostą i precyzyjną, bo mówi konkretnie czy firma jest w dobrej czy złej kondycji (w przeciwieństwie do analizy wskaźnikowej, która jest niejednoznaczna). Jej wadą jest jednak na stałe ustalenie wielkości parametrów (na podstawie testów na innych próbach z przeszłości). Podobną rolę do analizy dyskryminacyjnej odgrywają liniowy model prawdopodobieństwa oraz regresja logistyczna (model logitowy; logit), ale dostarczają one jeszcze dokładniejszej informacji: jakie jest prawdopodobieństwo, że na przykład firma jest w dobrej albo złej kondycji; jakie jest prawdopodobieństwo bankructwa. A zatem dostajemy nie tyle sygnał kupna lub sprzedaży, co prawdopodobieństwo sygnału (jakkolwiek to dziwnie brzmi). Należy dodać, że logit jest to metoda dużo bardziej skomplikowana i trudniejsza do zrozumienia i żeby ją zastosować konieczne jest użycie programu komputerowego, najlepiej do ekonometrii. Ja używam gretla.

Skoro chcemy uzyskać szansę tylko na to, że będzie dobrze lub źle, to znaczy, że zmienna Y ma tylko dwie wartości: 0 (będzie źle) lub 1 (będzie dobrze). Zacznijmy więc od najprostszego przykładu rzutu monetą, która nie musi być sprawiedliwa. Powiedzmy orzeł - 1, reszka - 0. Chcemy orła. Rzucamy raz monetą - udało się. Rzucamy drugi raz - znów orzeł. Za trzecim razem jest reszka. Rzucamy tak wielokrotnie i zapisujemy wyniki. Teraz chcemy sprawdzić prawdopodobieństwo wyrzucenia orła. Oczywiście nie potrzebujemy do tego żadnych cudów, logitów. Zwyczajnie obliczamy częstość wyrzuconych orłów i dzielimy przez liczbę wszystkich rzutów. Mamy więc to czego szukaliśmy. Powiedzmy, że zastępujemy monetę samochodem jadącym po trudnej drodze - 1 to brak wypadku, 0 to wypadek. Chcielibyśmy sprawdzić szansę, że wypadku nie będzie. Porównując auto z monetą napotykamy problem - nie możemy go testować wielokrotnie, ale tylko raz. Co zrobić? Moglibyśmy wprawdzie sprawdzić statystyki wypadków w kraju w ciągu roku i podzielić na liczbę jadących samochodów w tym samym czasie. Pomijając kwestię aktywnych samochodów, to taka metoda obliczania prawdopodobieństwa przypomina przysłowie o średniej liczbie nóg człowieka i psa. Przecież podejrzewamy, że aby doszło do wypadku, muszą zajść jakieś okoliczności, chcemy ustalić jakieś zmienne, które do tego prowadzą. Bez dokładniejszych informacji, obliczanie średniej częstości nie ma dla nas sensu. Chcemy większej precyzji, stąd pierwszy pomysł odpada. Musimy odkryć związek między jakimiś zmiennymi a zdarzeniem się wypadku, a następnie znaleźć jakiś sposób na obliczenie prawdopodobieństwa wypadku. W przypadku logitu będzie w zasadzie na odwrót: najpierw znajdziemy prawdopodobieństwo, a ono doprowadzi nas do związku między zmiennymi.

Liniowy model prawdopodobieństwa 

Skoro szukamy zmiennej zależnej i niezależnej, to powróćmy na chwilę do monety. Na wynik monety wpływała tylko jej budowa, a nie to, jak nią rzucamy. Budowa mogła powodować większe szanse wyrzucenia orła od reszki lub na odwrót. Co więcej, gdyby moneta była idealnie sprawiedliwa, to nawet budowa nie wpływałaby na zdarzenie oprócz tego, że pozwala jedynie na 0 lub 1. Innymi słowy w tym przypadku Y = X, gdzie Y to wyrzucenie 0 lub 1, a X to zmienna wpływająca na Y i przyjmuje to samo co Y. Jeżeli podejrzewamy, że moneta jest niesprawiedliwa, to to równanie musi się zmienić, np. na Y = a + b*X. Teraz moneta została jakby rozbita na dwie części: same zdarzenia Y oraz zdarzenia, które X wpływają na Y, ale nie są już tożsame z Y.
 
Wróćmy do przykładu z samochodem. Podobnie jak z monetą, musimy spytać co będzie miało wpływ na wypadek. Gdyby Y = X, to znaczyłoby, że na wypadek wpływa jedynie sama budowa samochodu i jego "wewnętrzne zagrożenia" albo sam kierowca. Głównym czynnikiem tworzącym zagrożenie są umiejętności kierowcy i jego stan psychiczny. Dla uproszczenia połączymy obie te cechy w jedną zmienną X: kierowca o wysokich umiejętnościach i w dobrym stanie psychicznym (1) lub nie spełniający tych wymagań (0). Gdyby nie było innych zakłócających czynników, to równanie Y = X mogłoby zostać spełnione: dobry kierowca (X = 1) nigdy nie miałby wypadku (Y = 1), natomiast zły kierowca (X = 0) zawsze miałby wypadek (Y = 0). Ta idealna korelacja rzecz jasna nie jest prawdziwa, ale pozwala dobrze zobrazować początkową sytuację. Znalezienie prawdopodobieństwa sprowadzałoby się do określenia z jakim kierowcą mamy do czynienia. I nagle spostrzegamy, że szansa musiałaby wynosić dokładnie 1 albo 0, bo utożsamiliśmy idealnie zdolności kierowcy z zagrożeniem wypadku. Skąd ta różnica w porównaniu z monetą? Rzucanie sprawiedliwej monety (bo o takiej teraz mówimy) to czysty los, natomiast zdolności i stan kierowcy losowe nie są. Można te zdolności oczywiście traktować jak zmienną losową, jednakże to właśnie fakt, że nie są w pełni losowe, pozwala na lepsze przewidywanie zdarzeń.

Oczywiście taka idealna zależność nie występuje, dlatego równanie Y = X należy także tutaj poprawić. Pytanie jak? Zacznijmy od liniowej zależności Y = a + b*X. Jeżeli X = 0, to Y = a, zaś gdy X = 1, to Y = a + b. To jednak ciągle oznacza doskonałą korelację równą 1: X to zawsze przekształcony Y. Załóżmy, że a = 0. Wtedy zawsze Y = 0 gdy X = 0. Dla b = 1, gdy X = 1, Y = 1. Zmienne zachowują się identycznie. W rzeczywistości relacja X-Y jest statystyczna, więc powinna być modelem regresji. To znaczy, że powinna zawierać dodatkowo składnik losowy. Czyli moglibyśmy zapisać

gdzie ε to składnik losowy. Wtedy dostajemy linię regresji jako warunkową wartość oczekiwaną Y względem X:

Zauważmy, że wartość oczekiwana Y równa się prawdopodobieństwu zajścia zdarzenia Y:

(1)

gdzie p to prawdopodobieństwo, że Y przyjmie wartość 1.

W tym miejscu ważne spostrzeżenie: skoro p to prawdopodobieństwo, że Y = 1, to znaczy, że:

(2)

 
Analogicznie warunkowa wartość oczekiwana Y będzie się równać prawdopodobieństwu warunkowemu zajścia zdarzenia Y pod warunkiem X:

(3)

Znów ważne jest to, że mamy tu do czynienia z prawdopodobieństwem warunkowym p, że Y wyniesie 1:

(4)

To kluczowy moment do zrozumienia dalszego wywodu. Warunkowa wartość oczekiwana Y to pewna średnia pod warunkiem, że X wyniesie (0 lub 1). A skoro średnia, to jest to tylko coś między 0 a 1. Ta średnia będzie jednak równa prawdopodobieństwu, że Y właśnie wyniesie 1 (pod warunkiem, że X równa się coś). To trudny, ale ważny element teorii. 

Jeżeli regresja Y byłaby liniowa, to dostajemy:

Taki model nazywa się liniowym modelem prawdopodobieństwa. Parametry a i b znajdujemy normalnie metodą MNK lub inną.

Przykład 1.
Stworzyłem przykładowy zestaw 50 danych: X - Czy kierowca ma zdolności? 1 - Tak, 0 - Nie. Y - Czy brak wypadku? 1 - Tak, 0 - Nie.

Wchodzimy do gretla. Model ->Klasyczna metoda najmniejszych kwadratów. Zaznaczamy jako zmienną niezależną Czy_kierowca_ma_zdolnosci, a jako zależną - Czy_bez_wypadku. Wynik:

 Z uzyskanego okna wyników wchodzimy w Prognozę. Otrzymujemy:

Model wskazuje teoretyczne prawdopodobieństwa (prognozy), jego błędy ex-ante i przedział ufności. Weźmy pierwszą obserwację: był wypadek (0), a prognoza, że nie będzie wypadku to 30%. Skąd wiem, że o tym mówi prognoza? Właśnie stąd, co wcześniej napisałem, że prognoza to prawdopodobieństwo, że Y = 1 pod warunkiem X. To pod warunkiem X powoduje, że nie powiedziałem wszystkiego. Prognoza, że nie będzie wypadku wynosi 30%, ale tylko dla tej jednej obserwacji. A ta jedna obserwacja spełnia pewne warunki. Tym warunkiem jest zmienna X. To znaczy, że dostajemy szansę 30%, że nie będzie wypadku, ale tylko wtedy, gdy X równa się... i tu musimy wrócić do pierwszego rysunku, gdzie jest lista wartości zmiennej Czy_kierowca_ma_zdolnosci. Pierwsza wartość wynosi zero. Czyli wiemy już, że gdy X = 0 (gdy kierowca nie ma zdolności) to prawdopodobieństwo, że Y = 1 (że nie będzie miał wypadku), wynosi tylko 0,3, co jest równoznaczne z tym, że będzie miał wypadek na 70% (1 - 0,3 = 0,7).

 Dla drugiej obserwacji Y = 1, a p = 0,6. Sprawdzamy pierwszy rysunek, wartość X(2) = 1. Czyli wiemy już, że gdy X = 1 (gdy kierowca ma zdolności) to prawdopodobieństwo, że Y = 1 (że nie będzie miał wypadku), wynosi 0,6.

Żeby pokazać zmienne niezależne i prognozę, zapisuję prognozy (+ w oknie prognozy) jako zmienną Pr_Czy_bez_wypadku_lin1. W oknie głównym zaznaczam tę nową zmienną razem z Czy_kierowca_ma_zdolnosci, i wyświetlam tabelę:

Faktycznie, dostaję to czego oczekiwaliśmy.

Podsumowując ten przykład, jeżeli mamy dobrego kierowcę, to spokojnie przejedziemy z prawdopodobieństwem 0,6, a jeżeli mamy złego kierowcę, to będziemy mieć wypadek z prawdopodobieństwem 0.7. Widzimy więc, że nie ma symetrii między odwrotnymi przypadkami.

Przykład 2.
Spróbujmy model poprawić, dodając drugi czynnik - otoczenie. Nowa zmienna niezależna X2 będzie pytać czy inni kierowcy mają zdolności? Tak - 1, Nie - 0. Dane tak samo własnoręcznie stworzyłem. Tabela wszystkich zmiennych z wartościami jest teraz taka:

Wszystko robimy tak samo jak poprzednio, tylko dodajemy nową zmienną niezależną do modelu.

Model się znacznie poprawił (AIC, BIC, HQC). Spójrzmy na prognozy:

Jak widać prognozy wypadków już dużo bardziej się zmieniają niż w pierwszym modelu. Mimo wszystko możemy bez nowych obserwacji sami wychwycić p. Zapiszmy prognozę w bazie danych (znak + w oknie prognozy), np. o nazwie Pr_Czy_bez_wypadku_lin2. W oknie głównym zaznaczmy zmienne niezależne i nowo dodaną zmienną. Wyświetlamy:

Gdy X1 = 0 i X2 = 0, to p = 0,13. Gdy X1 = 1 i X2 = 1, to p = 0,826. Gdy X1 = 1 i X2 = 0, p = 0,4. Gdy X1 = 0 i X2 = 1, p = 0,554. Te wielkości potem się powtarzają, stąd łatwo zgadniemy jakie będzie p dla nowych obserwacji (nowych kierowców). Jednak przy 3 zmiennych będzie to już na tyle skomplikowane, że nie opłaca się samemu tego robić i lepiej skorzystać z automatycznej prognozy.

Powiedzmy, że dostajemy nową obserwację X, tzn. nowego kierowcę, o którym wiemy, że bywa pijany, czyli X1 = 0, a do tego kierowcy, którzy jadą obok niego łamią prawo, bo stracili prawo jazdy, tj. X2 = 0. Nie wiemy, czy nastąpi wypadek, czyli nie posiadamy wartości Y. Dostajemy realną sytuację, w której przewidujemy przyszłość, tzn. chcemy się dowiedzieć, jaka jest szansa, że nie będzie wypadku. W Dane -> Dodaj obserwacje dodajemy 1 obserwację, zaznaczamy zmienną X1, X2 i Dane - > Edycja wartości. Możemy wpisać 51 obserwację, tj. X1 = 0, X2 = 0. Co ważne, nie dodajemy Y. Teraz ponownie wchodzimy w KMNK, zaznaczamy wszystko tak jak poprzednio i wchodzimy w prognozę. Dostałem:

Otrzymaliśmy to, czego oczekiwaliśmy, tj. p = 0,13. Jak wyżej napisałem, można samemu do tego dojść przy małej liczbie zmiennych, ale przy większej już taka procedura byłaby niezbędna. 
 
W ten sposób uzyskujemy prosty model do szacowania poziomu zagrożenia.

Liniowy model prawdopodobieństwa jest prosty i dość logiczny, ale ma dwie wady. Po pierwsze, w tablicy prognoz, przedział ufności zawiera liczby ujemne, a także powyżej 1. A przecież prawdopodobieństwo nie może być ujemne ani wyższe od 1. W omawianym przypadku problem ten ma charakter czysto teoretyczny, bo nie dostaniemy prognoz poniżej 0 ani powyżej 1. Niemniej mogą się pojawić przykłady, że gdy takie wartości się pojawią.  

Druga wada jest z punktu widzenia analizowanego przykładu poważniejsza. Bo na jakiej właściwie podstawie założyliśmy, że zależność Y-X jest liniowa? I stąd, dlaczego prawdopodobieństwo zdarzenia Y ma być funkcją liniową X? Gdy stosujemy klasyczną regresję liniową, to zazwyczaj patrzymy graficznie czy obserwacje X układają się liniowo względem Y. A jak się dzieje w tym przypadku? Wykres punktowy wygląda następująco:

Oczywiście wiele z 50 danych nakłada się na siebie, bo wszystkie mają tylko 0 lub 1.

Regresja logistyczna 

Czy można rozwiązać nasz problem bez założenia liniowej zależności X-Y? Powiedzmy, że mamy jakąś ogólną statystyczną zależność X od Y:

(5)

 gdzie ε to składnik losowy.

Rozwiązanie tego zadania wiąże się z wykorzystaniem zasady maksymalnej entropii. To ona doprowadza do wzoru na regresję logistyczną, który wyprowadziłem poprzednim razem:

(6)

gdzie:
p(i) - prawdopodobieństwo i-tego zdarzenia dla zmiennej X, która przyjmuje wartości x(i, 1), x(i, 2), ..., x(i, K).
Reszta to stałe.
Ponieważ są dwa zdarzenia Y, to i = 1, 2. Natomiast dla X jest w sumie K*i + 1 zdarzeń.

Mamy tu czynienia z K-wymiarową zmienną X. W naszym przypadku K = 1:

(7)

 W uproszczeniu:

(8)

Model (7) można też przekształcić do postaci:

(9)

Ponieważ szukamy p(Y = 1 | X), to (7), dla którego i = 1 (bo i = 1, 2), zapiszemy:

(10)

Pamiętajmy przy tym, że ciągle jest spełniona tożsamość (4), stąd:

(11)

Sumarycznie moglibyśmy zapisać:

Zmienną z formalnie należy traktować jako stałą wartość. Dla wszystkich zdarzeń, tj. X, p stanie się funkcją prawdopodobieństwa - dystrybuantą rozkładu logistycznego:

(12)

Powstaje nieco zagmatwana sprawa. Zgodnie z (12) zmienna Z ma rozkład logistyczny. Ale jedyną zmienną losową tworzącą Z jest X, zatem X też musi mieć rozkład logistyczny (tzn. rozkład prawdopodobieństwa zmiennej Z jest rozkładem logistycznym). Co się dzieje jednak z Y? Zmienna ta nie może być tożsama ze zmienną Z, bo by idealnie korelowała z X bez błędów losowych. To znaczy, że brakuje jej składnika losowego w Z. Pytanie czy ten składnik można po prostu dodać do Z? Okazuje się, że dla jednej zmiennej X można tak zrobić, ale nie gdy występuje więcej zmiennych. Np. jeśli X = 1, a beta(0) = 0, beta(1) = 0,5 to Y = 0,5 + składnik losowy. Wtedy E(Y | X) = 0,5, jeżeli składnik losowy jest niezależny od X i ma wartość oczekiwaną 0 (tak powinno być). Jeżeli jednak występują dwie zmienne niezależne, X1 i X2, wtedy jeśli beta(0) = 0, beta(1) = 0,5 oraz beta(2)  = 0,7, to E(Y | X) = 0,5 + 0,7 = 1,2, czyli p > 1. Do tego nie może dojść, gdyż nie dopuszcza do tego sama konstrukcja wzoru na p w (12). Prawdopodobieństwo po prostu nie może przekraczać 1. Wobec tego niemożliwe jest zwyczajne dodanie składnika losowego do Z, aby dostać Y. Oznacza to, że Y będzie funkcją nieliniową i poszukiwana funkcja g będzie postaci:

W sumie możemy jedynie napisać:

I tak jak stwierdziliśmy wyżej, w przypadku tylko jednej zmiennej niezależnej X - ale tylko zero-jedynkowej - prawdziwe będzie równanie:

co oznacza, że dla funkcji jednej zmiennej regresja logistyczna pokrywa się z modelem liniowym prawdopodobieństwa. Jaka jest tego przyczyna? Po prostu X ma tylko dwa punkty, a więc Y też musi mieć! Skacze więc z jednego punktu do drugiego i stąd sztucznie zachowuje się jak funkcja liniowa.
Dla zmiennej X nie-zerojedynkowej już tak nie będzie. 

Pozostało zadanie estymacji parametrów. Chociaż doszliśmy do dość zadziwiającego liniowego wyniku dla 1 zmiennej X, który prowadzi do wniosku, że można parametry poszukać MNK, to w ogólnym przypadku korzysta się z metody największej wiarygodności. Poszukuje się w niej w takich parametrów, aby prawdopodobieństwo obserwowanych danych (zmiennych o rozkładzie logistycznym) było największe. Na szczęście to trudne zadanie wykonuje program.


Przykład 3. 
Wykorzystamy te sam zestaw 50 danych co w przykładzie 1: X - Czy kierowca ma zdolności? 1 - Tak, 0 - Nie. Y - Czy brak wypadku? 1 - Tak, 0 - Nie:

Wchodzimy w gretla. Model -> Ograniczona zmienna zależna -> model logitowy -> binarny. Zaznaczamy jako zmienną niezależną Czy_kierowca_ma_zdolnosci, a jako zależną - Czy_bez_wypadku. Dodatkowo zaznaczamy opcję Pokaż wartość p (zamiast Pokaż efekt krańcowy dla średniej). Dostałem taki wynik:

Model okazuje się istotny statystycznie. Poniżej widzimy macierz błędów: Dane przewidywane (Y) kontra empiryczne (X): wiersz drugi, kolumna 2, czyli 1x1, oznacza ile razy model trafnie przewidywał brak wypadku (18). Mówiąc inaczej, ile razy występowało Y = 1 oraz X = 1. Wiersz pierwszy, kolumna pierwsza, czyli 0x0, oznacza ile razy model trafnie przewidywał wypadek (14). Mówiąc inaczej, ile razy występowało Y = 0 oraz X = 0. Pierwszy wiersz, druga kolumna, czyli 0x1 pokazuje ile razy model błędnie przewidywał, że kierowca przejedzie bez wypadku (12). Drugi wiersz, pierwsza kolumna, czyli 1x0 pokazuje ile razy model błędnie przewidywał, że kierowca będzie miał wypadek (6). Jak widać najwięcej błędów model popełnił przewidując brak wypadku, gdy kierowca miał zdolności, natomiast stosunkowo niedużo błędów popełnił szacując wystąpienie wypadku, gdy kierowca nie miał zdolności. Jest to sygnał, że słabe zdolności kierowcy rzeczywiście będą prognozować wypadek, natomiast wysokie zdolności mogą nie być wystarczającym czynnikiem do oceny czy nastąpi wypadek.

Zanim poprawimy model, sprawdźmy w końcu to, do czego od początku dążyliśmy - jakie jest prawdopodobieństwo Y = 1. Z uzyskanego okna wyników wchodzimy w Prognozę. Otrzymałem:

A zatem tak samo jak w liniowym modelu, każda obserwacja generuje prawdopodobieństwo wypadku, p(i, t). Jest to możliwe, bo wykorzystujemy teoretyczny poziom prawdopodobieństwa, dla którego, jak wyżej napisałem, X przyjmuje konkretne wartości dla t-tej obserwacji. Wg modelu, jeżeli kierowca ma zdolności, to szansa na brak wypadku, wynosi tylko 60%. A jeśli nie ma zdolności, to 30%. Czyli dostajemy to samo co w liniowym modelu, co wyżej sami przewidzieliśmy.

Przykład 4.
Spróbujmy model poprawić dodając nową zmienną niezależną X2, tę samą co w przykładzie 2: czy inni kierowcy mają zdolności? Tak - 1, Nie - 0. Tabela wszystkich zmiennych z wartościami jest teraz taka:

Wszystko robimy tak samo jak poprzednio, tylko dodajemy nową zmienną niezależną do modelu. Wynik dla logit:

Wszystkie kryteria, od BIC do AIC się poprawiły, czyli model się poprawił. Widać to zresztą po macierzy błędów: sumarycznie wcześniej było 18 błędów (12+6), teraz 14 (6+8).  Nowa zmienna okazuje się bardzo istotna stat. Teraz prognoza:

Zapiszmy prognozę w bazie danych (znak + w oknie prognozy), np. o nazwie Pr_Czy_bez_wypadku_log2. Wszystkie kolejne kroki możemy powtórzyć z przykładu 2. Czyli znów, chcąc prognozować wypadek, mając nowe informacje o X1 i X2, dopisujemy je do obserwacji i sprawdzamy jaka będzie prognoza Y.

Widzieliśmy, że dla X o 1 zmiennej prognozy w modelu liniowym i logicie były identyczne. Porównajmy teraz prognozy z przykładu 2 i przykładu 4:

Okazuje się, że prognozy są bardzo podobne, choć tym razem nieco się różnią.

Na koniec zaznaczę, że zarówno liniowy model prawdopodobieństwa jak i regresja logistyczna nie muszą być stosowane tylko dla zero-jedynkowej zmiennej X. To Y musi być zawsze binarny. Więcej o tym, z przykładami, napiszę następnym razem.