Regresja logistyczna w R: glm() z family = binomial

Gdy wynikiem jest tak/nie

Regresja liniowa przewiduje liczbę. Wiele pytań, które warto modelować, jest jednak binarnych: czy klient odejdzie, czy pacjent wyzdrowieje, czy ktoś kliknie w e-mail. Dopasowanie prostej do wyniku 0/1 od razu się sypie: prosta radośnie przewiduje prawdopodobieństwa −0.3 albo 1.4, co jest bezsensem.

Regresja logistyczna naprawia to, modelując prawdopodobieństwo wyniku przez transformację log-odds (logit): log(p / (1 − p)) = intercept + slope × x. Skala log-odds obejmuje całą oś liczbową, więc równanie liniowe pasuje tam naturalnie, a przekształcenie z powrotem ściska każde przewidywanie do przedziału (0, 1) wzdłuż znajomej krzywej w kształcie litery S. Cena tej sztuczki: współczynniki żyją w skali log-odds, a cała sztuka czytania regresji logistycznej polega na przekładaniu ich z powrotem na coś zrozumiałego dla ludzi.

Dopasowanie: glm() z family = binomial

glm() (uogólniony model liniowy) to starsze rodzeństwo lm(); family = binomial wybiera regresję logistyczną. W mtcars zmienna am zapisuje typ skrzyni biegów (1 = manualna, 0 = automatyczna). Czy oszczędne samochody częściej mają manualną skrzynię?

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

Dwie rzeczy przed odczytaniem wyniku. Po pierwsze, family = binomial nie jest opcjonalne: pomiń je, a glm() po cichu dopasuje zwykłą metodę najmniejszych kwadratów. Po drugie, zmienna wynikowa musi być binarna: 0/1, logiczna albo czynnik o dwóch poziomach (R modeluje prawdopodobieństwo drugiego poziomu).

A teraz summary, blok po bloku:

• Coefficients: Estimate dla mpg wynosi około 0.31 i jest to nachylenie w skali log-odds: każda dodatkowa mila na galon dodaje 0.31 do log-odds bycia manualem. Dodatni znak oznacza "zwiększa prawdopodobieństwo", ujemny "zmniejsza". Poza znakiem niczyja intuicja w tej skali nie działa i dlatego istnieje kolejna sekcja.
• z value i Pr(>|z|): ta sama logika co w testach t regresji (Estimate ÷ Std. Error, a potem p-value dla pytania "czy to może być zero?"), tylko z przybliżeniem normalnym, stąd z zamiast t. Tutaj p ≈ 0.011: związek między mpg a typem skrzyni raczej nie jest szumem.
• Null deviance kontra Residual deviance: dewiancja to w świecie glm miara złego dopasowania (mniej = lepiej). Null deviance (43.2 przy 31 df) dotyczy modelu z samym wyrazem wolnym, a residual deviance (29.7 przy 30 df) twojego modelu. Spadek o około 13.6 za cenę jednego stopnia swobody to w glm odpowiednik "R-kwadrat wzrósł".
• AIC: wynik do porównywania modeli, który równoważy dopasowanie i złożoność; niższy wygrywa. Sam w sobie nic nie znaczy, ale przydaje się między modelami kandydującymi na tych samych danych.

Od log-odds do ilorazów szans: exp(coef())

Podniesienie e do potęgi przenosi współczynniki z addytywnej skali log-odds do multiplikatywnej skali szans:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

exp(0.307) ≈ 1.36, a oto uczciwy wzorzec zdania do zapamiętania: "każda dodatkowa mila na galon mnoży szanse manualnej skrzyni biegów przez około 1.36." Iloraz szans powyżej 1 zwiększa szanse, poniżej 1 je zmniejsza, a dokładnie 1 oznacza brak efektu. Dlatego przy ilorazach szans werdykt z przedziału ufności brzmi "czy przedział wyklucza 1?" (a nie zero, bo zero było granicą w skali log-odds).

Uważaj na język: szanse to nie prawdopodobieństwa. Szanse = p / (1 − p), więc prawdopodobieństwo 0.75 to szanse 3. Pomnożenie szans przez 1.36 to nie to samo co pomnożenie prawdopodobieństwa przez 1.36, a przy częstym wyniku różnica jest duża. Iloraz szans 2 dla rzadkiego wyniku zachowuje się jak "mniej więcej podwojone ryzyko", ale dla wyniku występującego w 50% przypadków zdecydowanie nie. Nigdy nie opisuj ilorazu szans słowami ryzyka względnego ("1.36 razy bardziej prawdopodobne"), chyba że wynik jest rzadki.

Przewidywane prawdopodobieństwa: pułapka type = "response"

Najczęstszy błąd w regresji logistycznej spotykany w praktyce:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

Pierwsze wywołanie zwraca domyślne type = "link", czyli przewidywania w skali log-odds, razem z wartościami ujemnymi. Drugie zwraca prawdziwe prawdopodobieństwa. Jeśli twoje "prawdopodobieństwa" kiedykolwiek wychodzą ujemne albo większe niż 1, to właśnie z tego powodu. Uruchom blok: samochód z 15 mpg praktycznie nie ma szans na manualną skrzynię, samochód z 30 mpg bardzo prawdopodobnie ją ma, a krzywa S wygina się pośrodku.

Klasyfikacja: próg i tabela pomyłek

Prawdopodobieństwa stają się przewidywanymi klasami po wybraniu progu (domyślnie 0.5), a uczciwą kartą wyników jest tabela wartości przewidywanych i rzeczywistych:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

Komórki na przekątnej to trafne decyzje, a dwie komórki poza przekątną to dwa różne błędy (przewidywanie manualnej skrzyni dla automatu i odwrotnie). Sama ogólna trafność potrafi mocno schlebiać modelowi: jeśli 95% klientów nie odchodzi, "przewiduj, że nikt nie odejdzie" daje 95% i nie wyłapuje ani jednego odchodzącego. Dlatego zawsze patrz na oba typy błędów. A 0.5 to konwencja, nie prawo: gdy oba błędy mają różne koszty, przesuń próg odpowiednio.

Jedno uczciwe zastrzeżenie: ta tabela ocenia model na tych samych danych, na których go dopasowano, co mu schlebia. Prawdziwa ocena odkłada dane, których model nigdy nie widział.

Wiele predyktorów

Dokładnie jak w lm(): dodawaj składniki przez +, a każda interpretacja dostaje dopisek "przy pozostałych stałych":

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

Każdy współczynnik po podniesieniu e do potęgi to teraz mnożnik szans przy wzroście danego predyktora o jedną jednostkę wśród samochodów podobnych pod względem pozostałych predyktorów. Mechanizm się skaluje, ale skalują się też zastrzeżenia z regresji liniowej: skorelowane predyktory przetasowują sobie nawzajem współczynniki.

Ostrzeżenia

• Całkowita separacja. Jeśli predyktor idealnie rozdziela wynik (każdy samochód powyżej pewnego mpg jest manualny, każdy poniżej automatyczny), współczynnik największej wiarygodności dąży do nieskończoności. R ostrzega komunikatem glm.fit: fitted probabilities numerically 0 or 1 occurred i podaje ogromne współczynniki z absurdalnymi błędami standardowymi. Nie publikuj tych liczb: uprość model, zbierz więcej danych albo użyj metody z karą (pakiety brglm2 lub logistf).
• Wystarczająca liczba zdarzeń. Wiążącym ograniczeniem jest liczba rzadszych wyników, a nie łączna liczba wierszy. Stara reguła kciuka mówi o rzędzie 10-15 zdarzeń na predyktor. Przykłady z 32 samochodami służą tu do nauki mechaniki, a nie jako wzór wielkości próby do publikacji.
• Ilorazy szans to nie ryzyka względne, gdy wynik jest częsty. Było o tym wyżej, ale powtarzamy, bo recenzenci to wyłapią, nawet jeśli ty nie.

Co warto zapamiętać

• Wynik binarny → glm(y ~ x, data = df, family = binomial); nigdy nie zapominaj o family.
• Surowe współczynniki to log-odds; exp(coef(fit)) daje ilorazy szans, a wartością zerowego efektu dla ich przedziałów jest 1.
• Wzorzec zdania: "każdy wzrost x o jedną jednostkę mnoży szanse wyniku przez exp(b)."
• predict(..., type = "response") dla prawdopodobieństw: domyślne ustawienie zwraca log-odds, co jest pomyłką numer jeden.
• Klasyfikuj z progiem i oceniaj tabelą pomyłek; sama trafność potrafi kłamać.
• Uważaj na ostrzeżenia o separacji, licz swoje zdarzenia i nie przebieraj ilorazów szans za ryzyka względne.

Dalej: mechanizm stojący za każdym dotychczasowym przedziałem, czyli przedziały ufności z t.test(), confint() i prop.test().