Modele zmiennych jakościowych — logit i probit
Pełny wykład modeli dla binarnej zmiennej objaśnianej: wady liniowego modelu prawdopodobieństwa, wyprowadzenie logitu i probitu z modelu zmiennej ukrytej, log-szanse i iloraz szans z dowodem, estymacja największej wiarygodności i wklęsłość log-wiarygodności, efekty krańcowe — z wykresami w skali szarości (2D i 3D) oraz trzema pełnymi przykładami ekonomicznymi: zatrudnienie a wykształcenie, spłata kredytu a dochód, zakup a wiek — z kodem R i Python komentowanym linia po linii i bibliografią z DOI.
Wiele decyzji ekonomicznych ma charakter binarny: bezrobotny znajduje pracę albo nie, kredytobiorca spłaca dług albo bankrutuje, gospodarstwo kupuje dobro trwałe albo nie. Zmienna objaśniana przyjmuje wtedy wartości $\{0,1\}$, a celem jest modelowanie prawdopodobieństwa sukcesu. Artykuł przedstawia teorię logitu i probitu, a następnie stosuje ją do trzech realnych problemów ekonomicznych, przeprowadzając obliczenia dla każdego z nich kolejno.
Liniowy model prawdopodobieństwa i jego wady
Najprostszym podejściem jest zwykła MNK na zmiennej zero-jedynkowej — liniowy model prawdopodobieństwa (LPM): $P(Y_i=1\mid X_i)=\beta_0+\beta_1 X_i$. Jego zaletą jest prostota interpretacji ($\beta_1$ to zmiana prawdopodobieństwa), lecz ma trzy poważne wady.
Trzy wady LPM.
- Predykcje poza $[0,1]$ — dla skrajnych $X$ model przewiduje $P<0$ lub $P>1$, co jest bezsensowne.
- Wbudowana heteroskedastyczność — dla zmiennej zero-jedynkowej $\mathrm{Var}(Y\mid X)=P(1-P)$ zależy od $X$, więc założenie sferyczności jest z definicji złamane.
- Stały efekt krańcowy — LPM zakłada, że wpływ $X$ jest taki sam przy $P=0{,}05$ i $P=0{,}5$, co jest nierealistyczne; w rzeczywistości efekt powinien słabnąć przy krańcach.
Model zmiennej ukrytej i funkcje wiążące
Rozwiązaniem jest przepuszczenie kombinacji liniowej przez dystrybuantę, która utrzymuje wynik w $(0,1)$. Najczystsze uzasadnienie daje model zmiennej ukrytej.
Istnieje nieobserwowalna zmienna ukryta (użyteczność, skłonność) $y^\ast=X\beta+\varepsilon$, a obserwujemy jedynie jej znak:
$$ Y=\begin{cases}1,& y^\ast>0,\\ 0,& y^\ast\le0.\end{cases} $$Rozkład $\varepsilon$ wyznacza model: $\varepsilon$ logistyczne daje logit, $\varepsilon$ normalne — probit.
Przy symetrycznym rozkładzie błędu o dystrybuancie $F$:
$$ P(Y=1\mid X)=F(X\beta),\qquad F(z)=\begin{cases}\Lambda(z)=\dfrac{1}{1+e^{-z}}&\text{(logit)},\\[2mm]\Phi(z)&\text{(probit)}.\end{cases} $$- Definicja zdarzenia. $Y=1$ zachodzi dokładnie wtedy, gdy $y^\ast>0$, czyli $X\beta+\varepsilon>0$, równoważnie $\varepsilon>-X\beta$.
- Prawdopodobieństwo przez dystrybuantę. $P(Y=1\mid X)=P(\varepsilon>-X\beta)=1-F(-X\beta)$.
- Wykorzystanie symetrii. Dla rozkładu symetrycznego $1-F(-z)=F(z)$, więc $P(Y=1\mid X)=F(X\beta)$.
- Konkretne postacie. Dla $\varepsilon$ logistycznego $F=\Lambda$; dla $\varepsilon\sim\mathcal{N}(0,1)$ — $F=\Phi$. Obie funkcje są sigmoidalne, więc $P\in(0,1)$ przy malejącym efekcie krańcowym na krańcach.
Log-szanse i iloraz szans
W modelu logitowym logarytm szans jest liniowy w regresorach:
$$ \ln\frac{P(Y=1\mid X)}{P(Y=0\mid X)}=X\beta, $$a $e^{\beta_j}$ to iloraz szans: wzrost $x_j$ o jednostkę mnoży szanse sukcesu przez $e^{\beta_j}$.
- Szanse. Z $P=\Lambda(X\beta)$ mamy $1-P=\tfrac{e^{-X\beta}}{1+e^{-X\beta}}$, więc $\dfrac{P}{1-P}=e^{X\beta}$.
- Logarytm. $\ln\dfrac{P}{1-P}=X\beta$ — log-szanse są liniowe.
- Iloraz szans. Po wzroście $x_j$ o $1$ szanse zmieniają się z $O$ na $e^{\beta_j}O$, więc ich stosunek to $e^{\beta_j}$.
Estymacja największej wiarygodności
Dla niezależnych obserwacji log-wiarygodność wynosi
$$ \ell(\beta)=\sum_{i=1}^{n}\Big[y_i\ln F(X_i\beta)+(1-y_i)\ln\big(1-F(X_i\beta)\big)\Big]. $$Dla logitu równania wiarygodności mają postać $\sum_i\big(y_i-\Lambda(X_i\beta)\big)\mathbf{x}_i=\mathbf{0}$, a funkcja $\ell$ jest wklęsła, więc maksimum jest jedyne i wyznaczane numerycznie (Newton-Raphson).
- Pochodna funkcji logistycznej. $\Lambda'(z)=\Lambda(z)\big(1-\Lambda(z)\big)$.
- Gradient. Pochodna członu daje $(y_i-P_i)\mathbf{x}_i$ (z $P_i=\Lambda(X_i\beta)$); sumując, $\partial\ell/\partial\beta=\sum_i(y_i-P_i)\mathbf{x}_i$.
- Hesjan. $\partial^2\ell/\partial\beta\partial\beta^\top=-\sum_i P_i(1-P_i)\mathbf{x}_i\mathbf{x}_i^\top$.
- Wklęsłość. Każdy składnik jest dodatnio półokreślony ($P_i(1-P_i)>0$), więc Hesjan jest ujemnie półokreślony — $\ell$ jest wklęsła, a maksimum globalne.
Efekty krańcowe
Współczynnik $\beta_j$ nie jest efektem krańcowym, bo model jest nieliniowy. Efekt krańcowy zależy od punktu: $\dfrac{\partial P}{\partial x_j}=f(X\beta)\,\beta_j$ (dla logitu $f=\Lambda(1-\Lambda)$). Jest największy przy $P=0{,}5$ i maleje ku krańcom. Raportuje się zwykle AME (średni efekt krańcowy po całej próbie).
Przykład: wykształcenie a zatrudnienie
Kontekst i dane
Badamy, czy liczba lat nauki podnosi prawdopodobieństwo bycia zatrudnionym. Dla 12 osób obserwujemy wykształcenie (lata nauki) oraz status zatrudnienia (1 = pracuje, 0 = nie pracuje).
| Osoba | Lata nauki $x$ | Zatrudniony $y$ |
|---|---|---|
| 1 | 8 | 0 |
| 2 | 9 | 0 |
| 3 | 10 | 0 |
| 4 | 11 | 1 |
| 5 | 12 | 0 |
| 6 | 12 | 1 |
| 7 | 14 | 1 |
| 8 | 15 | 1 |
| 9 | 16 | 1 |
| 10 | 17 | 1 |
| 11 | 18 | 1 |
| 12 | 20 | 1 |
Krok 1 — Specyfikacja modelu
Ponieważ $y$ jest zero-jedynkowe, stosujemy logit: $P(\text{zatrudniony}=1\mid x)=\Lambda(\beta_0+\beta_1 x)$. Chcemy oszacować $\beta_0,\beta_1$ metodą największej wiarygodności.
Krok 2 — Estymacja (Newton-Raphson)
Startujemy z $\beta=(0,0)$ i iterujemy krok $\beta\leftarrow\beta-\mathbf{H}^{-1}\mathbf{g}$, gdzie $\mathbf{g}=\sum_i(y_i-P_i)\mathbf{x}_i$, a $\mathbf{H}=-\sum_i P_i(1-P_i)\mathbf{x}_i\mathbf{x}_i^\top$. Po kilku iteracjach uzyskujemy zbieżność:
$$ \hat\beta_0=-14{,}73,\qquad \hat\beta_1=1{,}303. $$Krok 3 — Iloraz szans
$$ e^{\hat\beta_1}=e^{1{,}303}=3{,}68. $$Każdy dodatkowy rok nauki mnoży szanse zatrudnienia przez $3{,}68$ (wzrost o $268\%$).
Krok 4 — Prawdopodobieństwa i próg
$$ P(\text{zatr.}\mid x{=}12)=\Lambda(-14{,}73+1{,}303\cdot12)=\Lambda(0{,}90)=0{,}71,\qquad P(\text{zatr.}\mid x{=}16)=\Lambda(6{,}11)=0{,}998. $$Próg $P=0{,}5$ wypada przy $x=-\hat\beta_0/\hat\beta_1=11{,}3$ roku nauki.
Krok 5 — Efekt krańcowy i dopasowanie
Średni efekt krańcowy $\mathrm{AME}=0{,}095$: przeciętnie każdy dodatkowy rok nauki podnosi prawdopodobieństwo zatrudnienia o około 9,5 punktu procentowego. Pseudo-$R^2$ McFaddena wynosi $0{,}64$ — bardzo dobre dopasowanie.
Wynik i interpretacja
Wykształcenie silnie i istotnie zwiększa szanse zatrudnienia. Osoba z 12 latami nauki ma około $71\%$ szans na pracę, a z 16 latami — praktycznie pewność ($99{,}8\%$). Granica, powyżej której zatrudnienie staje się bardziej prawdopodobne niż bezrobocie, to nieco ponad 11 lat nauki. Efekt krańcowy nie jest stały: dodatkowy rok nauki robi największą różnicę w okolicach progu (11–12 lat), a niewielką dla osób bardzo słabo lub bardzo dobrze wykształconych.
Kod R — pełne rozwiązanie
edu <- c(8,9,10,11,12,12,14,15,16,17,18,20) # lata nauki
zatr <- c(0,0,0,1,0,1,1,1,1,1,1,1) # 1 = zatrudniony
m <- glm(zatr ~ edu, family = binomial(link="logit")) # estymacja MNW (logit)
coef(m) # b0 = -14,73 ; b1 = 1,303
exp(coef(m)["edu"]) # iloraz szans = 3,68 na +1 rok nauki
predict(m, data.frame(edu=c(12,16)), type="response") # P = 0,71 ; 0,998
library(margins); summary(margins(m)) # AME = 0,095 (9,5 pkt proc.)
1 - m$deviance / m$null.deviance # pseudo-R2 McFaddena = 0,64
Kod Python — pełne rozwiązanie
import numpy as np, statsmodels.api as sm
edu = np.array([8,9,10,11,12,12,14,15,16,17,18,20.]) # lata nauki
zatr = np.array([0,0,0,1,0,1,1,1,1,1,1,1]) # 1 = zatrudniony
X = sm.add_constant(edu) # dodaj wyraz wolny
m = sm.Logit(zatr, X).fit() # estymacja MNW
print(m.params) # [-14.73, 1.303]
print(np.exp(m.params[1])) # iloraz szans = 3.68
print(m.predict([[1,12],[1,16]])) # P = 0.71 ; 0.998
print(m.get_margeff().summary()) # AME = 0.095
print(m.prsquared) # pseudo-R2 McFaddena = 0.64
Przykład: dochód a spłata kredytu
Kontekst i dane
Bank bada, czy miesięczny dochód kredytobiorcy wpływa na terminową spłatę kredytu. Dla 12 klientów mamy dochód (tys. zł/mies.) i wynik spłaty (1 = spłacony, 0 = niespłacony).
| Klient | Dochód $x$ [tys. zł] | Spłacony $y$ |
|---|---|---|
| 1 | 2,0 | 0 |
| 2 | 2,5 | 0 |
| 3 | 3,0 | 0 |
| 4 | 3,2 | 1 |
| 5 | 3,5 | 0 |
| 6 | 4,0 | 1 |
| 7 | 4,5 | 1 |
| 8 | 5,0 | 1 |
| 9 | 5,5 | 1 |
| 10 | 6,0 | 1 |
| 11 | 7,0 | 1 |
| 12 | 8,0 | 1 |
Krok 1 — Model
$P(\text{spłata}=1\mid x)=\Lambda(\beta_0+\beta_1 x)$, gdzie $x$ to dochód w tys. zł.
Krok 2 — Estymacja
Metodą największej wiarygodności (Newton-Raphson):
$$ \hat\beta_0=-12{,}61,\qquad \hat\beta_1=3{,}70. $$Krok 3 — Iloraz szans (kluczowa interpretacja)
$$ e^{\hat\beta_1}=e^{3{,}70}=40{,}5. $$Wzrost dochodu o 1 tys. zł mnoży szanse spłaty przez ponad 40. Tak wysoki iloraz szans wskazuje, że dochód jest zmienną o szczególnie silnej sile różnicującej w tej próbie.
Krok 4 — Prawdopodobieństwa
$$ P(\text{spłata}\mid x{=}4)=\Lambda(-12{,}61+3{,}70\cdot4)=\Lambda(2{,}19)=0{,}90. $$Próg $P=0{,}5$ wypada przy dochodzie $x=12{,}61/3{,}70=3{,}41$ tys. zł.
Krok 5 — Efekt krańcowy
$\mathrm{AME}=0{,}233$: przeciętnie każdy dodatkowy tysiąc złotych dochodu podnosi prawdopodobieństwo spłaty o około 23 punkty procentowe.
Wynik i interpretacja
Dochód jest bardzo silnym predyktorem spłacalności. Klient zarabiający 4 tys. zł ma $90\%$ szans terminowej spłaty; punkt równowagi (równe szanse spłaty i niespłaty) leży przy dochodzie około 3,4 tys. zł. Dla banku oznacza to, że dochód powyżej tego progu wyraźnie obniża ryzyko kredytowe — iloraz szans ($\times 40$ na każdy tysiąc) pokazuje, jak szybko ryzyko maleje wraz z dochodem. Uwaga: tak duży iloraz szans przy małej próbie ($n=12$) bywa niestabilny i wymaga ostrożności.
Kod R — pełne rozwiązanie
doch <- c(2,2.5,3,3.2,3.5,4,4.5,5,5.5,6,7,8) # dochod [tys. zl/mies.]
spl <- c(0,0,0,1,0,1,1,1,1,1,1,1) # 1 = kredyt splacony
m <- glm(spl ~ doch, family = binomial) # logit (link domyslny)
coef(m) # b0 = -12,61 ; b1 = 3,70
exp(coef(m)["doch"]) # iloraz szans = 40,5 na +1 tys. zl
predict(m, data.frame(doch=4), type="response") # P(splata) = 0,90
-coef(m)[1]/coef(m)[2] # prog P=0,5 przy dochodzie 3,41 tys. zl
library(margins); summary(margins(m)) # AME = 0,233
Kod Python — pełne rozwiązanie
import numpy as np, statsmodels.api as sm
doch = np.array([2,2.5,3,3.2,3.5,4,4.5,5,5.5,6,7,8.]) # dochod [tys. zl]
spl = np.array([0,0,0,1,0,1,1,1,1,1,1,1]) # 1 = splacony
m = sm.Logit(spl, sm.add_constant(doch)).fit() # estymacja MNW
print(m.params) # [-12.61, 3.70]
print(np.exp(m.params[1])) # iloraz szans = 40.5
print(m.predict([[1,4]])) # P(splata|4 tys.) = 0.90
print(-m.params[0]/m.params[1]) # prog P=0.5 -> 3.41 tys. zl
print(m.get_margeff().summary()) # AME = 0.233
Przykład 3 — Czy wiek wpływa na zakup produktu? (logit a probit)
Kontekst i dane
Sklep internetowy bada, czy wiek klienta wpływa na prawdopodobieństwo zakupu premium-produktu. Dla 12 klientów mamy wiek (lata) i decyzję (1 = kupił, 0 = nie kupił). Porównamy logit i probit.
| Klient | Wiek $x$ [lata] | Kupił $y$ |
|---|---|---|
| 1 | 20 | 0 |
| 2 | 25 | 0 |
| 3 | 28 | 1 |
| 4 | 30 | 0 |
| 5 | 33 | 1 |
| 6 | 35 | 1 |
| 7 | 38 | 0 |
| 8 | 40 | 1 |
| 9 | 45 | 1 |
| 10 | 50 | 1 |
| 11 | 55 | 1 |
| 12 | 60 | 1 |
Krok 1 — Dwa modele
Szacujemy logit $P=\Lambda(\beta_0+\beta_1 x)$ oraz probit $P=\Phi(\gamma_0+\gamma_1 x)$, by porównać wyniki.
Krok 2 — Estymacja logitu
$$ \hat\beta_0=-6{,}30,\qquad \hat\beta_1=0{,}203,\qquad e^{\hat\beta_1}=1{,}23. $$Krok 3 — Estymacja probitu
$$ \hat\gamma_0=-3{,}80,\qquad \hat\gamma_1=0{,}122. $$Krok 4 — Porównanie skali współczynników
Stosunek nachyleń $\hat\beta_1/\hat\gamma_1=0{,}203/0{,}122=1{,}67$ — bliski teoretycznemu $\approx1{,}6$. Mimo różnej skali współczynników przewidywane prawdopodobieństwa obu modeli są niemal identyczne.
Krok 5 — Efekt krańcowy
Dla logitu $\mathrm{AME}=0{,}026$: przeciętnie każdy dodatkowy rok wieku podnosi prawdopodobieństwo zakupu o około 2,6 punktu procentowego. Pseudo-$R^2$ McFaddena $=0{,}38$ — umiarkowane dopasowanie (wiek tłumaczy decyzję tylko częściowo).
Wynik i interpretacja
Starsi klienci kupują premium-produkt chętniej, lecz związek jest słabszy niż w poprzednich przykładach (umiarkowane pseudo-$R^2$): wiek to tylko jeden z wielu czynników decyzji. Iloraz szans $1{,}23$ oznacza, że każdy rok podnosi szanse zakupu o $23\%$, a efekt krańcowy ($\approx 2{,}6$ pkt proc. na rok) jest największy w wieku średnim, gdzie decyzja jest najbardziej niepewna. Praktyczny wniosek: wybór między logitem a probitem nie ma tu znaczenia merytorycznego — obie krzywe się pokrywają.
Kod R — pełne rozwiązanie
wiek <- c(20,25,28,30,33,35,38,40,45,50,55,60) # wiek [lata]
kup <- c(0,0,1,0,1,1,0,1,1,1,1,1) # 1 = kupil
mL <- glm(kup ~ wiek, family = binomial(link="logit")) # logit
mP <- glm(kup ~ wiek, family = binomial(link="probit")) # probit
coef(mL) # logit: -6,30 ; 0,203
coef(mP) # probit: -3,80 ; 0,122
coef(mL)["wiek"]/coef(mP)["wiek"] # stosunek nachylen = 1,67 (~1,6)
exp(coef(mL)["wiek"]) # iloraz szans = 1,23 na +1 rok
library(margins); summary(margins(mL)) # AME = 0,026 (2,6 pkt proc.)
Kod Python — pełne rozwiązanie
import numpy as np, statsmodels.api as sm
wiek = np.array([20,25,28,30,33,35,38,40,45,50,55,60.]) # wiek [lata]
kup = np.array([0,0,1,0,1,1,0,1,1,1,1,1]) # 1 = kupil
X = sm.add_constant(wiek)
mL = sm.Logit(kup, X).fit() # logit
mP = sm.Probit(kup, X).fit() # probit
print(mL.params) # [-6.30, 0.203]
print(mP.params) # [-3.80, 0.122]
print(mL.params[1]/mP.params[1])# 1.67 (~1.6)
print(np.exp(mL.params[1])) # iloraz szans = 1.23
print(mL.get_margeff().summary()) # AME = 0.026
Ocena dopasowania
Metoda MNW nie daje klasycznego $R^2$. Stosuje się: pseudo-$R^2$ McFaddena ($1-\ell/\ell_0$; wartości $0{,}2$–$0{,}4$ uchodzą za dobre), tablicę trafności (odsetek poprawnych klasyfikacji przy progu $0{,}5$), krzywą ROC i pole AUC oraz test ilorazu wiarygodności. Dla więcej niż dwóch kategorii stosuje się logit wielomianowy lub uporządkowany.
Wkład polskiej szkoły
Modele logitowe i probitowe znalazły w Polsce szerokie zastosowanie w prognozowaniu upadłości przedsiębiorstw — nurcie rozwijanym m.in. przez Dariusza Hadasika i Elżbietę Mączyńską, którzy konstruowali modele dyskryminacyjne i logitowe dla polskich danych finansowych. Kanonicznym opracowaniem metodyki jest Mikroekonometria pod redakcją Marka Gruszczyńskiego (SGH).
Podsumowanie
- LPM: prosty, lecz przewiduje poza $[0,1]$, jest heteroskedastyczny i ma stały efekt krańcowy.
- Zmienna ukryta: $Y=1\iff X\beta+\varepsilon>0$; $P=F(X\beta)$, logit ($\Lambda$) lub probit ($\Phi$).
- Log-szanse: $\ln\tfrac{P}{1-P}=X\beta$; iloraz szans $e^{\beta_j}$.
- MNW: log-wiarygodność wklęsła; Newton-Raphson.
- Efekty krańcowe: $f(X\beta)\beta_j$, zależne od punktu; raportuj AME.
- Logit ≈ probit: $\beta_{\text{logit}}\approx1{,}6\,\beta_{\text{probit}}$; predykcje niemal identyczne.
Bibliografia
- Berkson, J. (1944). Application of the Logistic Function to Bio-Assay. Journal of the American Statistical Association, 39(227), 357–365. https://doi.org/10.1080/01621459.1944.10500699
- Bliss, C. I. (1934). The Method of Probits. Science, 79(2037), 38–39. https://doi.org/10.1126/science.79.2037.38
- Cramer, J. S. (2003). Logit Models from Economics and Other Fields. Cambridge University Press.
- McFadden, D. (1974). Conditional Logit Analysis of Qualitative Choice Behavior. W: P. Zarembka (red.), Frontiers in Econometrics (s. 105–142). Academic Press.
- Gruszczyński, M. (red.). Mikroekonometria. Modele i metody analizy danych indywidualnych. Wolters Kluwer.
Dalej: Testy hipotez · Zmienne instrumentalne · Dane panelowe
- Podręcznik: J. M. Wooldridge, Introductory Econometrics: A Modern Approach, rozdz. 17
- J. S. Cramer, Logit Models from Economics and Other Fields, Cambridge University Press, 2003
- M. Gruszczyński (red.), Mikroekonometria. Modele i metody analizy danych indywidualnych, Wolters Kluwer
- Pakiety: R
glm(family=binomial),margins; Pythonstatsmodels(Logit, Probit)