Lab 7 - regresja liniowa cz. 1

Sportowcy

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

• Ile ma obserwacji, a ile zmiennych? Które zmienne są ilościowe, a które jakościowe?
• Jaka jest średnia oraz wariancja każdej ze zmiennych ilościowych?
• Które zmienne ilościowe są najbardziej skorelowane, które najsłabiej, a które mają najsilniejszą korelację ujemną?

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()

Jak widać, zależność jest w przybliżeniu liniowa.
Korzystając ze wzorów podanych na początku tego skryptu, obliczamy parametry $\beta$.
W naszym przypadku $Y$ to waga Wt, a $X$ to wzrost Ht.

beta1 <- sum((Ht-mean(Ht))*(Wt - mean(Wt)))/sum((Ht - mean(Ht))^2)
beta0 <- mean(Wt) - beta1*mean(Ht)

Możemy przedstawić wynik na wykresie:

ggplot(ais, aes(x=Ht, y=Wt)) + geom_point() + geom_abline(slope=beta1, intercept=beta0, col='red')

Interpretacja współczynników jest następująca: Jeśli porównamy sportowców różniących się wzrostem o 1 cm, to średnia różnica ich wagi będzie równa $= $ 1.1171169. Sportowiec o wzroście 180 cm średnio waży $\beta_0 + 180\beta_1$ = 74.8920324 kg. Jeśli natomiast weźmiemy sportowca o zerowym wzroście, to będzie on średnio ważył $\hat{\beta_0} =$ -126.189011 kilo.

Ostatnie zdanie jest oczywiście bez sensu.
Pomimo tego, że parametr $\beta_0$ oznacza dokładnie to, co napisałem, to na ogół nie interpretujemy go w ten sposób.
Powodem jest to, że badane zależności są liniowe tylko w przybliżeniu i w ograniczonym przedziale. Gdybyśmy uwzględnili na powyższym wykresie dzieci, to zależność wagi od wzrostu przestałaby być liniowa.
Kosztem takiego przybliżenia jest to, że niektóre parametry tracą swój pierwotny sens. Parametr $\beta_0$ jest w tym przypadku sztucznym parametrem pozwalającym na lepsze dopasowanie się do danych. Dlatego, aby poprawnie interpretować otrzymane wyniki, trzeba dobrze rozumieć techniki, których się używa.

Zajmijmy się teraz statystyczną istotnością parametru $\hat{\beta_1}$. Chcemy zweryfikować hipotezę $H_0 : \beta_1 = 0$. Jeśli uda nam się ją odrzucić, to będziemy mogli stwierdzić, że dane wspierają założenie o istnieniu zależności pomiędzy wzrostem a wagą. W przeciwnym wypadku będziemy musieli uznać, że dane są zbyt słabe żeby stwierdzić jakąkolwiek zależność, że otrzymaliśmy niezerową wartość $\beta_1$ całkowicie przypadkiem, i że równie dobrze waga może (choć nie musi) być totalnie niezależna od wzrostu. Krótko mówiąc, że dane są do kitu.

Żeby przeprowadzić test, musimy założyć, że estymator $\hat{\beta_1}$ jest zmienną losową o pewnym rozkładzie prawdopodobieństwa (uwaga, estymator $\hat{\beta_1}$, a nie parametr $\beta_1$!).
Na wykładzie i ćwiczeniach zostało pokazane, że jeśli $\varepsilon \sim \mathcal{N}(0, \sigma^2)$, to $\beta_i \sim \mathcal{N}(\beta_i, \sigma^2 v_i)$, gdzie $v_i$ to element diagonali macierzy $(X^TX)^{-1}$ odpowiadający parametrowi $\beta_i$. Uwaga! Aby zachować zgodność z numeracją współczynników $\beta_0, \beta_1, \dots, \beta_k$, elementy $v_i$ również numerujemy od zera.

Zakładając, że $\hat{\beta}$ ma rozkład normalny, możemy przeprowadzić test t Studenta wykorzystując statystykę $$t = \frac{\hat{\beta_1}}{\hat{\sigma}\sqrt{v_1}},$$

gdzie $\hat{\sigma}^2 = \frac{1}{n - k - 1} \sum_{i=1}^n (Y_i - \hat{Y}_i)^2$ jest estymatorem nieobciążonym wariancji błędu $\varepsilon$, a $\hat{Y}_i$ to przybliżenie i-tej obserwacji $Y$ za pomocą funkcji liniowej (odpowiadające czerwonej linii na wykresie powyżej). Obserwowane wartości błędu losowego $\varepsilon$, czyli $Y_i - \hat{Y}_i$, nazywają się residuami.

Wykorzystamy teraz powyższe wzory, aby sprawdzić statystyczną istotność naszego parametru $\beta_1$. W pierwszym kroku obliczymy residua. Na ich podstawie ocenimy, czy możemy przyjąć, że $\varepsilon$ ma rozkład normalny.

predykcja_Wt <- beta0 + beta1 * Ht 
r <- Wt - predykcja_Wt
qqnorm(r)

Wykres kwantylowy w większości układa się na prostej, również po zbliżeniu na wybrane fragmenty (por. przykład zadłużenia gmin). Jest kilka nieco odstających punktów, ale nie powinno to zbyt mocno zaszkodzić naszym testom. Na wszelki wypadek można te punkty usunąć - decyzję pozostawiam Tobie.

Teraz możemy wyestymować wariancję błędu. Zwróć uwagę, że estymator nieobciążony wariancji błędu ma mniejszy mianownik w porównaniu z estymatorem nieobciążonym wariancji ze zwykłej próbki statystycznej: mamy $n - k - 1$ zamiast spotykanego dotychczas $n - 1$. Intuicyjny powód jest następujący: ponieważ prosta $\hat{Y} = X \hat{\beta}$, którą przeprowadzamy przez nasze punkty, minimalizuje błąd kwadratowy, to prawdziwa prosta $Y = X\beta$ będzie miała na ogół nieco większy błąd. Kwadraty residuów $(Y - X \hat{\beta})^2$ zaniżają zatem prawdziwe błędy kwadratowe $(Y - X\beta)^2 = \varepsilon^2$, dlatego zwiększamy je odpowiednio modyfikując mianownik.

sigma2 <- sum(r^2)/(length(Ht)-2)  # u nas n = dlugosc wektora Ht, k = 1
sigma <- sqrt(sigma2)

Odchylenie standardowe błędu losowego ma naturalną interpretację: Znając wzrost sportowca, możemy przewidzieć jego wzrost jako $\beta_0 + \beta_1 X$, i wówczas średnio rzecz biorąc pomylimy się o 8.720304 kg. Równoważnie: po uwzględnieniu wzrostu, pozostało nam do wyjaśnienia jeszcze 8.720304 kg wagi. Jeśli odchylenie jest równe zero, to znaczy że przeprowadziliśmy prostą dokładnie przez wszystkie punkty, a więc wyjaśniliśmy całą informację zawartą w zmiennej $Y$.

Obliczymy teraz statystykę testową. Potrzebujemy do tego i-tej wartości na przekątnej macierzy $(X^TX)^{-1}$. Ponieważ w naszym przypadku macierz $X^TX$ jest wymiaru 2x2, to moglibyśmy zrobić to ręcznie. Dobrze jednak wiedzieć, w jaki sposób zrobić to w przypadku ogólnym, ponieważ ilustruje to bardzo częsty problem w matematyce stosowanej. Odwracanie macierzy jest numerycznie bardzo trudne i złożone obliczeniowo, dlatego zawsze należy tego unikać. Zamiast tego, należy starać się przekształcić problem odwracania macierzy do problemu rozwiązywania układu równań liniowych.

Niech $e_i$ będzie pionowym wektorem długości $n$ z jedynką na i-tej współrzędnej i zerami poza nią. Wówczas iloczyn $e_i^T (X^TX)^{-1} e_i$ jest równy i-temu wyrazowi na przekątnej macierzy $(X^TX)^{-1}$ (mnożenie przez $e_i^T$ z lewej daje nam i-ty wiersz, a przez $e_i$ z prawej i-tą kolumnę macierzy). Chcemy zatem obliczyć $v_i = e_i^T (X^TX)^{-1} e_i$. Oznaczmy $w = (X^TX)^{-1} e_i$. Wektor $w$ możemy teraz obliczyć, rozwiązując układ równań liniowych $X^TX w = e_i$. Następnie, aby znaleźć $v_i$, wystarczy wybrać i-tą współrzędną wektora $w$.

Zobaczmy jak to zaimplementować w praktyce aby obliczyć wartość $v_1$ (pamiętaj, że w tym przypadku wartości na diagonali numerujemy od 0). Na początku stwórzmy macierz planu $X$, zawierającą kolumnę złożoną z jedynek oraz kolumnę zawierającą zmienną Ht (jeśli nie pamiętasz, skąd wzięła się tu kolumna z jedynkami, to przypomnij sobie slajd 13 z wykładu 6). Następnie obliczymy macierz $X^TX$, korzystając z funkcji t() zwracającej macierz transponowaną oraz operatora %*%, który wykonuje mnożenie macierzowe.

X <- cbind(1, Ht)
XTX <- t(X) %*% X  

Teraz obliczymy wektor $w$, korzystając z funkcji solve():

e_1 <- c(0, 1)
w <- solve(XTX, e_1)

Ostatecznie bierzemy drugą współrzędną wektora $w$, równą 5.2502156\times 10^{-5} (pierwsza współrzędna odpowiada $v_0$). Możesz sprawdzić, że jest ona równa wartości obliczonej poprzez ręczne odwrócenie macierzy, czyli $v_1 = (\sum X_i^2 - n \bar{X}^2 )^{-1}$.

Na końcu obliczamy statystykę $t$ i obliczamy p-wartość, zakładając, że $t$ ma rozkład t Studenta z $n-k-1$ stopniami swobody.

statystyka_t <- beta1/(sigma*sqrt(w[2]))
p_wartosc <- 2*pt(-statystyka_t, length(Ht)-2) 

Otrzymujemy p-wartość równą 9.6386472\times 10^{-43}, wskazującą, że nasze dane bardzo silnie świadczą o istnieniu zależności pomiędzy wzrostem a wagą.

Skonstruujemy teraz przedział ufności dla parametru $\beta_1$. Z wykładu wiemy, że ma on postać

$$(\hat{\beta_1} - z_{1 - \alpha/2} \sqrt{v_1} \hat{\sigma}, \hat{\beta_1} + z_{1 - \alpha/2} \sqrt{v_1} \hat{\sigma}).$$

Korzystamy z obliczonych wcześniej wartości:

alpha <- 0.05
z <- qnorm(1 - alpha/2)
c <- z * sqrt(w[2]) * sigma
przedzial_ufnosci <- c(beta1 - c, beta1 + c)

Przedział ufności na poziomie 0.95 dla średniej zależności pomiędzy wzrostem w centymetrach a wagą w kilogramach jest równy (0.9932748, 1.240959). Ponieważ przedział ten jest stosunkowo wąski, to możemy mieć dość dużą pewność co do wartości tego parametru.

Po zakończeniu obliczeń należy jeszcze wywołać komendę detach(ais), która działa odwrotnie do funkcji attach() i schowa kolumny tablicy ais, co pomaga utrzymać porządek w kodzie.

detach(ais)

Ręczne przeprowadzenie takich obliczeń jak powyżej jest bardzo kształcące i warto to zrobić raz w życiu.
Ale w praktyce oczywiście wszystko robimy korzystając z gotowych funkcji w R.
Do stworzenia modelu liniowego służy funkcja lm() (skrót od linear model).
Funkcja ta jako pierwszy argument przyjmuje nowy typ zmiennej, czyli formułę. Jest to bardzo wygodny sposób określania modelu: napis Y ~ X oznacza, że chcemy modelować zmienną objaśnianą Y za pomocą zmiennej objaśniającej X. Podobnie jak w bibliotece ggplot2, nie używamy apostrofów wokół nazw zmiennych.
Jako drugi argument funkcji lm() podajemy ramkę danych z których funkcja weźmie zmienne do modelu.
Aby stworzyć model liniowy wyjaśniający wagę za pomocą wzrostu wystarczy zatem napisać:

model <- lm(Wt ~ Ht, ais)

Najważniejsze informacje o naszym modelu możemy wyświetlić za pomocą funkcji summary():

summary(model)

## 
## Call:
## lm(formula = Wt ~ Ht, data = ais)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.372  -5.296  -1.197   4.378  38.030 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -126.18901   11.39656  -11.07   <2e-16 ***
## Ht             1.11712    0.06319   17.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.72 on 200 degrees of freedom
## Multiple R-squared:  0.6098, Adjusted R-squared:  0.6079 
## F-statistic: 312.6 on 1 and 200 DF,  p-value: < 2.2e-16

W polu ‘Residuals’ mamy podsumowanie rozkładu residuów $Y - \hat{Y}$ (to te same residua, które obliczyliśmy wcześniej ręcznie jako r <- Wt - predykcja_Wt).
W polu ‘Coefficients’ mamy podsumowanie współczynników $\hat{\beta}$. Intercept to wyraz wolny, czyli $\hat{\beta_0}$. Kolumna ‘Estimate’ zawiera wartości estymatorów, kolumna ‘Std. Error’ ich odchylenia standardowe $\hat{\sigma}\sqrt{v_i}$, kolumna ‘t value’ wartości statystyki $t_i = \hat{\beta_i}/(\hat{\sigma}\sqrt{v_i})$, a kolumna ‘Pr(>|t|)’ zawiera prawdopodobieństwo otrzymania większej wartości statystyki niż wartość zaobserwowana, czyli po prostu p-wartości. Trzy gwiazdki za p-wartością oznaczają “bardzo istotny” wynik, w tym sensie, że mamy duże podstawy aby przypuszczać, że dany parametr nie jest równy zero.
Na końcu podsumowania mamy podane kilka przydatnych statystyk, do których jeszcze wrócimy.

Aby otrzymać przedział ufności dla parametrów $\beta$, wystarczy teraz użyć funkcji confint():

confint(model)

##                    2.5 %      97.5 %
## (Intercept) -148.6618436 -103.716178
## Ht             0.9925209    1.241713

Jednym z największych zalet regresji liniowej w porównaniu z wieloma innymi technikami modelowania jest możliwość otrzymania przedziałów ufności dla predykcji. Możemy nie tylko przewidzieć wagę sportowca o wzroście 180 cm, ale również ocenić, jak daleko od wartości przewidywanej może potencjalnie być wartość prawdziwa. W tym celu korzystamy z funkcji predict(), której podajemy model oraz ramkę danych z wartościami zmiennych dla którch chcemy wykonać predykcję. Parametr interval='prediction' sprawia, że funkcja zwróci przedział ufności dla predykcji na poziomie 95%, czyli taki przedział, który z prawdopodobieństwem 0.95 zawiera prawdziwą wartość zmiennej Wt odpowiadającą Ht=180.

predict(model, data.frame('Ht'=180), interval='prediction')

##        fit      lwr      upr
## 1 74.89203 57.65398 92.13008

Widzimy, że sportowcy o wzroście 180 kg ważą na ogół od 57.6 do 92.13 kg, przy czym średnio 74.89 kg.
Korzystając z funkcji predict(), możemy łatwo narysować przedział ufności dla predykcji na wykresie.

przedzial_ufnosci <- predict(model, data.frame('Ht'=seq(120, 230, by=1)), interval='prediction')
przedzial_ufnosci <- as.data.frame(przedzial_ufnosci)
przedzial_ufnosci$Wt <- seq(120, 230, by=1)

ggplot(ais, aes(x=Ht, y=Wt)) + geom_point() + geom_abline(slope=beta1, intercept=beta0, col='red') + geom_ribbon(aes(x=Wt, ymin=lwr, ymax=upr), data=przedzial_ufnosci, alpha=0.1)

Zadanie 3. Regresja wieloraka. Za pomocą regresji liniowej zbadaj zależność wagi sportowca od wszystkich pozostałych zmiennych ilościowych. Które zmienne mają najsilniejszy wpływ na procent tkanki tłuszczowej? Czy lepiej ocenić to na podstawie kolumny ‘Estimate’, ‘Std. Error’, ‘t value’, czy p-wartości?

Jak duży wpływ na wagę mają czynniki, których nie obserwujemy w danych (np. czynniki genetyczne)?
Jaka jest średnia różnica pomiędzy wagą dwóch sportowców, którzy różnią się wzrostem o 1 cm, a wszystkie pozostałe cechy mają identyczne? Czy ta różnica jest większa, czy mniejsza niż w poprzednim zadaniu?

Wskazówka. Jeśli chcemy modelować zmienną $Y$ jako kombinację liniową zmiennych $X1$ i $X2$, to piszemy formułę Y ~ X1 + X2. Jeśli chcemy wykorzystać wszyskie zmienne, to piszemy Y ~ . (kropka oznacza “wszystkie kolumny tabeli poza Y”). Jeśli chcemy wykluczyć zmienną $Z$, możemy napisać Y ~ . - Z.

Zmienne jakościowe

Za pomocą regresji liniowej można również badać, w jaki sposób na zmienną objaśnianą wpływają cechy jakościowe, kodowane jako factory.
Dla każego poziomu factora poza poziomem bazowym (zwanym również poziomem referencyjnym) tworzona jest osobna zmienna zero-jedynkowa, która przyjmuje wartość 1 gdy dana obserwacja odpowiada temu poziomowi.
Na przykład, jeśli chcemy uwzględnić w naszym modelu płeć, to w macierzy planu $X$ tworzymy dodatkową kolumnę, która przyjmie wartość 1 gdy dany sportowiec jest mężczyzną i 0 gdy jest kobietą (ponieważ poziom referencyjny kolumny ais$Sex to female). Jeśli używamy funkcji lm(), to zmienne jakościowe piszemy w formułach tak samo, jak ilościowe. Inny jest natomiast ich opis w podsumowaniu modelu - znajdzie się tam średnia wartość zmiennej $Y$ dla każdego z poziomów factora.

Zadanie 4. Stwórz model liniowy wyjaśniający wagę za pomocą wszystkich zmiennych poza zmienną Sport. Zinterpretuj wyniki. Czy z modelu wynika, że sportowcy różnych płci różnią się wagą? Dlaczego otrzymaliśmy taki rezultat?
Czy taki model przewiduje, że zależność pomiędzy wzrostem a wagą może być różna dla różnych płci?

Wskazówka. Funkcja summary() wyświetla informacje dotyczące zmiennych jakościowych w polu ‘Coefficients’ osobno dla każdego poziomu. Nazwa odpowiedniego współczynnika powstaje przez sklejenie nazwy zmiennej i nazwy poziomu. W naszym przypadku w polu ‘Coefficients’ znajduje się wiersz ‘Sexmale’, odpowiadający poziomowi ‘male’ kolumny ‘Sex’ z danych ‘ais’. Pole ‘Estimate’ mówi, w jaki sposób zmieni się waga sportowca, jeśli stanie się on mężczyzną przy zachowaniu wszystkich pozostałych cech na ustalonym poziomie.

Kwartet Anscombe’a

W ocenie modelu liniowego pomoże nam statystyka $R^2$, zwana współczynnikiem determinacji, i zdefiniowana jako wariancja zmiennej $Y$ wyjaśniona przez nasz model: $$R^2 = \frac{\sum_{i=1}^n (\hat{y}_i - \bar{y}_i)^2}{\sum_{i=1}^n (y_i - \bar{y}_i)^2}.$$ W podsumowaniu modelu liniowego, otrzymanym za pomocą funkcji summary(), możemy znaleźć wartość statystyki $R^2$ w polu ‘Multiple R-squared’.

Jak zobaczymy w następnym zadaniu, statystyka $R^2$ służy głównie do wykrywania słabych modeli. Wysoka wartość tej statystyki nie jest wystarczającym argumentem za tym, że model dobrze wyjaśnia dane, mimo tego, że wiele osób tak ją stosuje.

Zadanie 5. Załaduj dane anscombe za pomocą komendy data(anscombe). Dane mają postać macierzy z kolumnami $x_i$ oraz $y_i$ dla i = 1, 2, 3, 4. Następnie:

• Utwórz cztery modele liniowe y_i ~ x_i.
• Porównaj współczynniki determinacji tych modeli. Współczynniki możesz znaleźć korzystając z funkcji summary.
  
• Przedstaw zależność pomiędzy $y_i$ a $x_i$ na wykresie, i zaznacz na nim prostą regresji. Czy wszystkie modele są tak samo dobre?