Na dzisiejszych zajęciach przypomnimy sobie, jakie są główne założenia regresji liniowej, nauczymy się wykrywać odstępstwa od tych założeń, oraz zobaczymy, jak sobie z nimi radzić.

Diagnostyka modelu liniowego

Regresja liniowa opiera się na kilku założeniach, które łatwo zapamiętać jako zasadę LINE:

Regresja liniowa natomiast nie zakłada, że każdej wartości \(Y\) odpowiada tylko jedna wartość \(X\), ani na odwrót. Jest wręcz pożądane, żeby dla jednej wartości zmiennej objaśniającej \(X\) otrzymać kilka niezależnych pomiarów zmiennej objaśnianej \(Y\) (na przykład zmierzyć kilkukrotnie tempo reakcji chemicznej w tej samej temperaturze), ponieważ daje nam to więcej informacji o błędach losowych.
Jest to bardzo częste źródło błędów u studentów (i nie tylko), i spotkałem się wielokrotnie z uśrednianiem wartości \(Y\) odpowiadających tym samym wartościom \(X\). Nie ma żadnej potrzeby, żeby tak robić, a wręcz może to prowadzić do błędnych wyników.

Oczywiście, tak samo jak w przypadku testów statystycznych, w praktyce nic nie jest idealne i poza nielicznymi sytuacjami żadne z tych założeń tak naprawdę nie jest spełnione. Im bliżej jesteśmy do ich spełnienia, tym bardziej wiarygodne wyniki otrzymamy.

Odstępstwa od poszczególnych założeń zasady LINE mają bardzo różny wpływ na nasze wnioski.
W bardzo zgrubny sposób podsumowuje je poniższa tabela. Przy pierwszym czytaniu będzie prawdopodobnie bardzo niejasna, ale mam nadzieję, że po przerobieniu zadań będzie lepiej.

Założenie Konsekwencja odstępstwa Co można z tym zrobić
Liniowość trendu To podstawowe założenie całej metody, więc kiedy nie jest spełnione, to jest dosyć słabo. Parametry \(\beta\) tracą w zasadzie jakiekolwiek znaczenie, tzn. przestają oznaczać średnią zmianę \(Y\) gdy zmienimy \(X\) o jednostkę. Na szczęście istnieje dużo metod przekształcenia danych, które często pozwalają uzyskać liniowy trend.
Niezależność błędów Wszystko zależy od tego, w jaki sposób błędy są zależne. Może się popsuć wszystko albo nic. Wszystko zależy od struktury błędów. Dla pewnych specyficznych przypadków istnieją specjalne techniki przekształcania danych, na przykład analizowanie różnic pomiędzy kolejnymi obserwacjami.
Normalność błędów Dopóki wartość średnia błędów wynosi 0, to estymacja \(\beta\) jest poprawna. Nie działają testy statystyczne ani przedziały ufności. Czasami udaje się odpowiednio przekształcić dane, na przykład transformacją Boxa-Coxa, ale dość często po prostu trzeba się z tym pogodzić.
Homoskedastyczność Obszary na których jest wysoka wariancja błędów mają mocniejszy wpływ na estymację \(\beta\), ponieważ dużo silniej wpływają na odległość kwadratową \(Y\) od \(\hat{Y}\). Możemy przez to mieć większy błąd estymacji. Nie działają testy przedziały ufności. Wszystko zależy od struktury wariancji. Bardzo często wariancja rośnie wraz ze wzrostem \(X\) - wówczas może się udać przekształcić dane tak, aby uzyskać homoskedastyczność.

Jak widać, przekształcanie danych jest bardzo ważnym (i często trudnym) elementem regresji liniowej.
Możemy transformować zarówno zmienną objaśnianą \(Y\), jak i zmienne objaśniające \(X\).
Wbrew pozorom, w regresji liniowej daje to różne wyniki: modelowanie \(\log(Y) \sim X\) to nie to samo, co modelowanie \(Y \sim \exp(X)\)!

Rozpatrzmy pewną odwracalną funkcję \(f\) (na ogół \(f(x) = log(x)\) lub \(f(x) = x^\alpha\)), i załóżmy, że:

Nie ma żadnego gotowego przepisu w jaki sposób przekształcić dane, i często trzeba posługiwać się intuicją lub wiedzą dziedzinową. Zilustrujemy teraz o co w tym wszystkim chodzi na prostym przykładzie.

Zadanie 1. Zadanie przykładowe. Załaduj dane bats.csv, znajdujące się na stronie przedmiotu. Dane zawierają długości ramion nietoperzy oraz ich wiek w momencie pomiaru. Zbadaj zależność pomiędzy wiekiem nietoperza a długością jego ramienia. W jakim wieku ramię nietoperza jest rozwinięte w 90%?

Rozwiązanie. W pierwszym kroku, jak zawsze, wczytujemy dane i wizualizujemy je używając biblioteki ggplot2. Zmienna \(Y\) będzie oznaczać długość ramienia w milimetrach, a \(X\) wiek nietoperza w dniach.

library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
bats <- read.csv('bats.csv', header=T)
ggplot(bats, aes(x=days, y=forearm.length)) + geom_point() + theme_minimal()

Widzimy, że zależność nie jest liniowa. Wygląda za to z grubsza jak \(Y = \log(X)\). Zobaczmy zatem, co się stanie, jeśli przestawimy \(Y\) w zależności od \(\log(X)\):

ggplot(bats, aes(x=log(days), y=forearm.length)) + geom_point() + theme_minimal()

Zależność po transformacji wygląda na bardziej liniową niż wcześniej. Transformacja logarytmiczna ma dodatkowo tę zaletę, że jest łatwo interpretowalna: długość ramienia zwiększy się o 1 mm, gdy wiek wzrośnie około 2.7 raza, ponieważ \(log(X\codot e) = log(X) + 1 = Y + 1\).

Inną popularną transformacją jest transformacja Boxa-Coxa, w której dla pewnego parametru \(\lambda\) robimy przekształcenie \[ Y_i \rightarrow \frac{Y_i^\lambda - 1}{\lambda} \] lub \[ X_i \rightarrow \frac{X_i^\lambda - 1}{\lambda}. \] Dla \(\lambda \rightarrow 0\) dostajemy przekształcenie logarytmiczne, tzn \(Y_i \rightarrow \log(Y_i)\). Parametr \(\lambda\) jest dobrany tak, aby rozkład residuów był jak najbliższy rozkładowi normalnemu.

Funkcja boxcox z pakietu MASS wykonuje transformację Boxa-Coxa zmiennej zależnej \(Y\). Jej składnia jest bardzo podobna do funkcji lm():

library(MASS)
bc <- boxcox(forearm.length~days, data=bats, lambda=seq(0, 4, by=1/10))

Powyższy wykres przedstawia wartość logarytmu funkcji wiarygodności residuów (która ocenia podobieństwo do rozkładu normalnego) w zależności od parametru \(\lambda\). Optymalny parametr wybieramy następująco:

best_lambda <- bc$x[which.max(bc$y)]

Następnie przekształcamy naszą zmienną i wizualizujemy wynik:

bats$forearm_boxcox <- (bats$forearm.length**best_lambda - 1)/best_lambda
ggplot(bats, aes(x=days, y = forearm_boxcox)) + geom_point() + theme_minimal()

Wygląda na to, że transformacja Boxa-Coxa na zmiennej objaśnianej nie dała w tym przypadku zbyt dobrych rezultatów. Oczywiście, podobne podejście można stosować również do transformacji logarytmicznej, i badać log-wiarygodność na przykład dla przekształcenia postaci \(X \rightarrow \log(X + \lambda)\). Wymaga to jednak większego nakładu pracy. Niejednokrotnie nad jednym modelem liniowym można spędzić kilka tygodni.

Żeby tego uniknąć, na potrzeby tych zajęć zbudujemy model liniowy w oparciu o transformację logarytmiczną, a następnie przeprowadzimy diagnostykę modelu, aby ocenić, na ile odstępujemy od założeń LINE.

bats$log_days <- log(bats$days)
model <- lm(forearm.length ~ log_days, data=bats)
summary(model)
## 
## Call:
## lm(formula = forearm.length ~ log_days, data = bats)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1069 -1.8870  0.2385  2.0376  5.0147 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.7853     0.8988   10.89 6.12e-13 ***
## log_days      6.7094     0.3791   17.70  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.745 on 36 degrees of freedom
## Multiple R-squared:  0.8969, Adjusted R-squared:  0.8941 
## F-statistic: 313.2 on 1 and 36 DF,  p-value: < 2.2e-16

Współczynnik \(R^2\) jest bardzo wysoki, więc model wyjaśnia dużo wariancji w danych. Na podstawie współczynników możemy uznać, że \(Y \approx 9.78 + 6.70\log(X)\). To oznacza, że średnia długość ramienia jednodniowego nietoperza (\(X=1\)) wynosi około 9.78 mm, następnie rośnie o 6.70 mm gdy jego wiek zwiększa się \(e\)-krotnie (\(e \approx 2.72\)). Zwizualizujmy nasz model:

predykcja <- data.frame('days'=bats$days, 'forearm.length' = 9.78 + 6.70*log(bats$days))
ggplot() + geom_point(aes(x=days, y=forearm.length), data=bats) + geom_line(aes(x=days, y=forearm.length), data=predykcja, col='red') + theme_minimal()

Przyjrzymy się teraz wykresom diagnostycznym, które otrzymamy wywołując komendę plot(model). Po wpisaniu tej komendy, wciśnięcie klawisza Enter wyświetli kolejne wykresy. Możemy też skorzystać z argumentu which, który kontroluje, który wykres należy wyświetlić. Więcej informacji możemy uzyskać w dokumentacji funkcji plot.lm, którą uzyskamy wywołując komendę ?plot.lm.

plot(model, which=1)

Powyższy wykres przedstawia zależność \(Y\) od \(\hat{Y}\), czyli przybliżenia \(Y\) funkcją liniową (czyli predykcji \(Y\)). Służy on do oceny, czy zależność w danych jest liniowa. Na podstawie wykresu można uznać, że zależność jest w przybliżeniu liniowa, chociaż w praktyce prawdopodobnie szukalibyśmy lepszej transformacji.

Ten wykres pozwala również pobieżnie ocenić czy błędy są niezależne. Częstym źródłem zależności pomiędzy błędami jest autokorelacja - czyli korelacja \(\varepsilon_{i+1}\) z \(\varepsilon_{i}\) lub wcześniejszymi błędami. Na ogół interesuje nas korelacja dodatnia. Sprawia ona, że różnica pomiędzy dwoma sąsiadującymi residuami jest mniejsza niż pomiędzy dwoma wybranymi losowo. Takie zjawisko widać zresztą na powyższym wykresie. W tym przypadku jest ono jednak spowodowane tym, że modelowana przez nas zależność nie jest do końca liniowa. W zadaniach dodatkowych znajduje się konkretny przykład danych w których występuje typowa autokorelacja błędów.

plot(model, which=2)