Na dzisiejszych laboratoriach będziemy ćwiczyć regresję liniową, czyli technikę pozwalającą na badanie zależności liniowych pomiędzy zmiennymi losowymi. Jest to jedna z najważniejszych technik matematyki stosowanej, używana w niemal wszystkich innych dziedzinach nauki.

Regresja liniowa została omówiona na Wykładzie 6. Na początku podamy kilka informacji uzupełniających wiedzę z wykładu. Jeśli na początku ten paragraf wyda się niejasny, to warto do niego wrócić ponownie po przerobieniu zadań z dzisiejszego labu.

Jak wiadomo, ogólny wzór na współczynniki regresji liniowej to \(\hat{\beta} = (X^TX)^{-1}X^TY\), gdzie \(X\) to macierz planu zawierająca zmienne objaśniające (predyktory), a \(Y\) to wektor zawierający zmienną objaśnianą.
W przypadku jednej zmiennej objaśniającej, czyli gdy \(X\) jest wektorem zawierającym współrzędne \(X_1, \dots, X_n\), powyższy wzór na parametry \(\beta\) upraszcza się następująco: \[ \hat{\beta_1} = \frac{\sum_{i=1}^n (X_i - \bar{X})(Y_i - \bar{Y})}{\sum_{i=1}^n (X_i - \bar{X})^2} = \text{Cor}(X, Y) \cdot \text{sd}(Y) / \text{sd}(X) \] \[\hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{X}\] Zwróć uwagę, że we wzorze na \(\hat{\beta_1}\) możemy zidentyfikować trzy komponenty: Korelację, mierzącą stopień zależności pomiędzy \(X\) a \(Y\); Odchylenie standardowe \(Y\), mierzące, jak bardzo zmienia się ta zmienna; oraz odchylenie standardowe zmiennej \(X\). Zgodnie z intuicją, współczynnik \(\hat{\beta_1}\) jest tym większy, im większy jest rozrzut zmiennej \(Y\), i im mniejszy jest rozrzut zmiennej \(X\).

W przypadku wielu zmiennych objaśniających, wyrażamy zmienną \(Y\) jako kombinację liniową zmiennych objaśnianych oraz tak zwanego błędu losowego \(\varepsilon\): \(Y = \beta_0 + \beta_1 X_1 + \dots + \beta_k X_k + \varepsilon\). Nazwa błąd, choć często stosowana, jest tutaj dość niefortunna. Zmienna \(\varepsilon\) reprezentuje “efekty losowe” wpływające na zmienną \(Y\), ale w tym przypadku “losowość” oznacza tylko to, że nie jesteśmy w stanie wyjaśnić tych efektów za pomocą zmiennych \(X_1, \dots, X_k\). Innymi słowy, nasze zmienne objaśniające pozwalają nam wyjaśnić część zmienności zmiennej \(Y\), a zmienna \(\varepsilon\) oznacza zmienność której nie udało nam się wyjaśnić. Błąd losowy nie oznacza natomiast ani tego, że \(Y\) jest w jakiś sposób “niepoprawna”, ani tego, że jest “nieprzewidywalna” - wzięcie dodatkowych zmiennych objaśniających może zredukować błąd losowy do zera.

Wspomniane powyżej czynniki, których nie uwzględniliśmy w macierzy planu, modelujemy matematycznie jako zmienną losową \(\varepsilon\), ponieważ pozwala nam to na uzyskanie dodatkowych informacji o zależności pomiędzy zmiennymi \(Y\) i \(X\).
Rozkład zmiennej \(\varepsilon\) mówi nam o tym, jakie wartości może przyjmować zmienna \(Y\), jeśli zmierzymy ją przy ustalonych wartościach zmiennych \(X_1, \dots, X_k\) (dla przykładu z wykładu: jeśli sprawdzimy liczbę sprzedanych jednostek produktu w nowym sklepie o znanych nakładach na reklamę).

Losowość z \(\varepsilon\) “przenosi się” na zmienną \(Y\), więc tę ostatnią również traktujemy jako zmienną losową. Mamy m.in. \(\mathbb{E}(Y) = \beta_0 + \beta_1 X_1 + \dots + \beta_k X_k\) (przy założeniu że \(\varepsilon\) ma średnią 0).

Przy interpretacji modelu liniowego trzeba dobrze zdawać sobie sprawę z tego, co oznaczają współczynniki \(\beta\). Zmienna \(\beta_k\) jest równa średniej zmianie zmiennej \(Y\), jeśli zwiększymy zmienną \(X_k\) o 1, trzymając wszystkie pozostałe zmienne objaśniające na ustalonym poziomie: \(\beta_0 + \beta_1 X_1 + \dots + \beta_k (X_k + 1) = \mathbb{E}(Y) + \beta_k\). Wartość współczynnika mierzy zatem, jak silnie dana zmienna objaśniająca wpływa średnio na zmienną objaśnianą. Oczywiście, to nie oznacza, że \(Y\) zawsze zmieni się dokładnie o \(\beta_k\) jeśli zwiększymy \(X_k\) o 1, ponieważ grają tu rolę jeszcze czynniki losowe \(\varepsilon\).

Więcej na temat regresji liniowej można przeczytać m.in. na tej stronie. Przejdziemy teraz do zadań, które zilustrują powyższy paragraf. W pierwszym zadaniu wczytamy i wstępnie zbadamy zbiór danych na którym będziemy pracować.

Sportowcy

Zadanie 1. Wczytaj i obejrzyj zbiór danych dotyczący sportowców, który znajdziesz tutaj. Odpowiedz na poniższe pytania:

Utwórz wykresy punktowe obrazujące zależność między każdą parą zmiennych ilościowych. Możesz wykorzystać do tego funkcję pairs lub (bardziej estetyczną) funkcję ggpairs z pakietu GGally. Pokoloruj punkty wg płci. Jeśli wykres jest mało czytelny, wybierz podzbiór zmiennych ilościowych.

Poniższy wykres, otrzymany za pomocą funkcji ggpairs z biblioteki GGally, przedstawia zalezności pomiędzy czterema zmiennymi: Grubością fałdów skóry SSF, procentem tkanki tłuszczowej w ciele X.Bfat, ilością czerwonych krwinek Hc, i wagą Wt.
Wykresy na przekątnej obrazują rozkłady poszczególnych zmiennych. Górna część wykresu zawiera korelacje pomiędzy parami zmiennych, a dolna przedstawia zależności pomiędzy nimi za pomocą wykresu punktowego. Punkty są pokolorowane w zależności od płci.
Czy widzisz coś dziwnego w zależności pomiędzy wagą a poziomem tkanki tłuszczowej? A pomiędzy wagą i ilością czerwonych krwinek?

ggpairs(ais, aes(col=Sex), columns=c(9, 10, 5, 13))  

Zadanie 2 Zadanie przykładowe. Wykorzystaj regresję liniową, aby zbadać zależność wagi sportowców (Wt) od ich wzrostu (Ht). Sprawdź, czy zależność jest statystycznie istotna i znajdź przedział ufności na poziomie 95% dla współczynnika \(\beta_1\). Wykorzystaj wzory podane na slajdach do Wykładu 6. Następnie porównaj swoje wyniki z otrzymanymi za pomocą funkcji lm(). Wykorzystaj funkcję predict(), aby otrzymać przedział ufności na poziomie 95% dla sportowca o wzroście 180 cm.

Rozwiązanie.

W pierwszym kroku wczytamy dane. Ponieważ interesują nas wyłącznie kolumny Ht oraz Wt, to warto utworzyć zmienne które zapewnią bezpośredni dostęp do tych kolumn. Dzięki temu nie będziemy musieli za każdym razem pisać ais$Ht.

W R są na to dwa sposoby. Jeden z nich to przypisanie kolumny na nową zmienną x <- ais$Ht. Drugim sposobem jest wykorzystanie funkcji attach(), która sprawi, że kolumny ramki danych będą widoczne jako zmienne. Możemy dzięki temu obliczyć średnią kolumny Ht pisząc po prostu mean(Ht).

ais <- read.table("ais.txt", header = T)
attach(ais)
mean(Ht)
## [1] 180.104

Pierwszym krokiem analizy danych powinno być zawsze zwizualizowanie danych.
Zwizualizujmy zatem:

library(ggplot2)
ggplot(ais, aes(x=Ht, y=Wt)) + geom_point()