Regresja liniowa w R: lm(), summary() i predict()

Idea: prosta najmniejszych kwadratów

Regresja liniowa dopasowuje prostą do chmury punktów: y = intercept + slope × x. Spośród wszystkich możliwych prostych lm() wybiera tę, która minimalizuje sumę kwadratów reszt, gdzie reszta to pionowa odległość między punktem a prostą. Podnoszenie do kwadratu sprawia, że duże chybienia liczą się nieproporcjonalnie mocno, i dlatego jedna dzika obserwacja odstająca potrafi przechylić całe dopasowanie.

Korelacja daje jedną liczbę bez jednostek na pytanie "jak ściśle te zmienne poruszają się razem", a regresja daje równanie: z jednostkami, nachyleniem, które da się zinterpretować, i mechanizmem do przewidywania.

fit <- lm(mpg ~ wt, data = mtcars)
print(coef(fit))

Czytaj formułę jako "modeluj mpg jako funkcję wt". Dwa współczynniki to dopasowana prosta: mpg ≈ 37.3 − 5.3 × masa. Nachylenie ma prawdziwe jednostki: każde dodatkowe 1000 funtów masy samochodu (wt jest w jednostkach po 1000 funtów) kosztuje około 5.3 mili na galon. Wyraz wolny (37.3 mpg przy masie zero) to tylko miejsce, w którym prosta przecina zero. Żaden samochód nie waży zero, więc nie przywiązuj do niego zbyt dużej wagi.

Omówienie summary()

summary(fit) to wynik, który każdy kurs statystyki każe interpretować. Uruchom go, a potem przejdź blok po bloku:

fit <- lm(mpg ~ wt, data = mtcars)
summary(fit)

Call: powtarza dopasowany model. Teraz to drobiazg, ale ratuje skórę, gdy żonglujesz sześcioma obiektami modeli.

Residuals: pięcioliczbowe podsumowanie tego, co zostało (wartość rzeczywista − przewidywana). Mediana powinna być blisko 0, a Min/Max oraz 1Q/3Q mniej więcej symetryczne. Wyraźna asymetria sugeruje, że modelowi z prostą czegoś brakuje.

Coefficients: serce wyniku, jeden wiersz na każdy składnik:

• Estimate: dopasowana wartość. Dla wt to −5.34: każde dodatkowe 1000 funtów wiąże się z mniej więcej 5.3 mpg mniej. Zawsze przekładaj nachylenie na zdanie z jednostkami; to zdanie jest całą praktyczną treścią modelu.
• Std. Error: jak bardzo oszacowanie chwiałoby się przy powtarzanych próbach. Oszacowania oddalone od zera o najwyżej kilka błędów standardowych są niepewne.
• t value: Estimate ÷ Std. Error, czyli o ile błędów standardowych współczynnik jest oddalony od zera (tutaj −9.56).
• Pr(>|t|): p-value dla pytania "czy ten współczynnik może naprawdę wynosić zero?". Dla wt wynosi około 1.3e-10: gdyby masa naprawdę nie miała liniowego związku z mpg, tak strome nachylenie praktycznie nigdy nie pojawiłoby się w próbie 32 samochodów. Małe p-value to dowód, że związek istnieje, ale nie dowód, że model jest poprawny, i nie miara ważności (maleńki, precyzyjnie oszacowany efekt też dostaje maleńkie p-value).
• Signif. codes / gwiazdki: wizualny skrót kolumny z p-value. Wygodny, ale nie dodaje żadnej informacji.

Residual standard error: 3.05 on 30 degrees of freedom: typowy rozmiar błędu przewidywania w jednostkach zmiennej objaśnianej. Przewidywania mylą się zwykle o około 3 mpg. Oceniaj to na tle skali mpg (która mieści się mniej więcej w przedziale 10-34).

Multiple R-squared: 0.75: masa wyjaśnia około 75% wariancji mpg. Adjusted R-squared (0.74) przelicza to z karą za każdy predyktor, bo surowa wersja może tylko rosnąć, gdy dodajesz zmienne, nawet losowy szum. Przy porównywaniu modeli z różną liczbą predyktorów uczciwa jest wersja skorygowana. Oprzyj się też odruchowi "dobry model = wysokie R²": naprawdę przydatny efekt może żyć w modelu z niskim R² (zaszumiony wynik, jeden z wielu czynników), a wysokie R² może wynikać z przeuczenia albo z wycieku zmiennej.

F-statistic: 91.4 ... p-value: 1.29e-10: test całego modelu, czyli czy ten model jest lepszy niż "przewiduj wszystkim średnią". Przy jednym predyktorze powtarza test t dla nachylenia (zauważ, że 9.56² ≈ 91.4); przy kilku predyktorach staje się łącznym testem, że co najmniej jeden współczynnik jest niezerowy. Jego mechanizm to ten sam rozkład wariancji co w ANOVA.

Regresja wieloraka: przy pozostałych stałych

Predyktory dodajesz przez +:

fit2 <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit2)

Interpretacja zmienia się w jednym kluczowym punkcie. Każde Estimate to teraz efekt danego predyktora przy pozostałych predyktorach stałych: współczynnik wt (około −3.9, wcześniej −5.3) to koszt dodatkowej masy w mpg przy porównaniu samochodów o tej samej mocy. Wartość −5.3 z prostej regresji po cichu zawierała w sobie to, że cięższe samochody zwykle są też mocniejsze; regresja wieloraka to rozdziela. Z tego samego powodu współczynniki się przesuwają, gdy dodajesz zmienne: jeśli nowy predyktor koreluje ze starym, zmienia się zakres obowiązków starego. R-kwadrat rośnie do około 0.83, a tutaj uczciwym porównaniem z modelem z jednym predyktorem jest wersja skorygowana.

Przewidywanie: predict()

Dopasowany model to funkcja, a predict() ją oblicza. Zbuduj ramkę danych newdata, której nazwy kolumn dokładnie pasują do predyktorów:

fit <- lm(mpg ~ wt, data = mtcars)

new_cars <- data.frame(wt = c(2.5, 3.5))
print(predict(fit, newdata = new_cars))

print(predict(fit, newdata = new_cars, interval = "confidence"))
print(predict(fit, newdata = new_cars, interval = "prediction"))

Dwa typy przedziałów odpowiadają na różne pytania, a mylenie ich to klasyczny błąd na egzaminie:

• interval = "confidence": niepewność co do średniej, czyli "dla wszystkich samochodów ważących 2500 funtów, gdzie leży średnie mpg?". Wąski i zwęża się wraz z przyrostem danych.
• interval = "prediction": zakres, w którym prawdopodobnie znajdzie się pojedynczy nowy samochód o tej masie. Dużo szerszy, bo pojedynczy samochód ma własny rozrzut wokół prostej, którego nie uśredni żadna ilość danych.

Podawanie przedziału ufności, gdy pytanie dotyczy jednej nowej obserwacji, drastycznie zawyża twoją precyzję.

Diagnostyka i pułapka ekstrapolacji

summary() mówi, co model szacuje, a wykresy reszt mówią, czy mu wierzyć. W sesji interaktywnej:

par(mfrow = c(2, 2))
plot(fit)   # four diagnostic plots

Na co patrzeć: Residuals vs Fitted powinien być bezkształtną chmurą. Krzywa oznacza, że zależność nie jest prostoliniowa; lejek (rozrzut rosnący z wartościami dopasowanymi) oznacza niestałą wariancję, a wtedy twoje błędy standardowe są przekłamane. Punkty na wykresie Q-Q powinny przylegać do prostej: ciężkie ogony oznaczają, że obserwacje odstające zniekształcają dopasowanie. Scale-Location to ponownie test lejka. Residuals vs Leverage wskazuje punkty wpływowe, czyli obserwacje, które same ciągną współczynniki (w mtcars pojawiają się tu zwykle egzotyczne auta, np. Chrysler Imperial). Szybki wykres punktowy surowych danych przed dopasowaniem wyłapuje większość tych problemów wcześnie.

Na koniec pułapka, której nie wyłapie żadna diagnostyka: ekstrapolacja. Model uczył się na samochodach ważących mniej więcej 1500-5400 funtów. Podaj predict() wartość wt równą 8, a radośnie zwróci ujemne mpg: matematyka przedłuża prostą w nieskończoność, ale dowody kończą się na krawędzi danych. Przewiduj tylko w zakresie, na którym dopasowano model (albo blisko niego).

Co warto zapamiętać

• fit <- lm(y ~ x, data = df) dopasowuje prostą najmniejszych kwadratów; coef(fit) to równanie, summary(fit) pełny raport.
• Czytaj Estimate jako zdania z jednostkami; Pr(>|t|) pyta "czy to może być zero?", a nie "czy to ma znaczenie?".
• Residual standard error to typowe chybienie w prawdziwych jednostkach; adjusted R-squared to uczciwa liczba do porównywania modeli.
• W regresji wielorakiej każdy współczynnik oznacza "przy pozostałych stałych", a współczynniki przesuwają się, gdy dochodzą skorelowane predyktory.
• predict(fit, newdata, interval = ...): "confidence" dla średniej, "prediction" dla jednego nowego przypadku (ten szeroki).
• Sprawdzaj plot(fit) pod kątem krzywych, lejków i punktów wpływowych; nigdy nie ufaj przewidywaniom poza zakresem danych.

Dalej: gdy wynikiem jest tak/nie zamiast liczby, czyli regresja logistyczna z glm().