---
title: "Wielokrotna (wieloraka) regresja liniowa"
output:
  pdf_document: default
  html_document:
    df_print: paged
editor_options:
  chunk_output_type: console
---

Na dzisiejszych zajęciach opowiemy sobię trochę o wielokrotnej (wielorakiej) regresji liniowej (and. *multiple linear regression*). Pracować będziemy na przykładzie z rozdziału 15 podręcznika Howella, który znajduje się na stronie w materiałach dostępnych do pobrania po zalogowaniu. Proszę Państwa o przeczytanie tego rozdziału z podręcznika. Kod, który wklejam poniżej, ma stanowić bazę do ćwiczeń i ilustracje, jak różne rzeczy omówione w książce możemy zrobić w R. Kod będzie objaśniony, ale nie zastąpi (!) to lektury podręcznika. Dodatkowo pomyślałem, że możnaby przemycić trochę ciekawych i przydatnych R-owych bibliotek oraz technik, więc będę czasami korzystał z nieznanych Państwu funkcji, ale proszę się nie martwić -- wszystkie będę objaśniał i tłumaczył co robią.

Na początek wczytajmy dwie biblioteki, których użyjemy do stworzenia tego pliku oraz dane znajdujące się w pliku `Tab15-1.dat`. 

Biblioteka `pander` służy do ładnego wyświetlania wyników różnych testów statystycznych. Bez niej *output* różnych operacji, które wykonujemy, wygląda jak wydruk kodu źródłowego programu (szarawe tło oraz czcionka o stałej szerokości). Pander sprawia, że zamiast wydruku testy statystyczne zwracają ładnie sformatowane wyniki w formie np. tabelki. Z pakietu `knitr` będziemy używać funkcji `kable`, która w zasadzie służy do tego samego, czyli do ładnego wyświetlania danych w dokumentach HTML (takich jak ten). Od razu mogą Państwo zobaczyć działanie tego pakietu -- dzięki funkcji `kable` nasza ramka danych jest estetycznie sformatowana.

Uwaga! Wczytaliśmy dane za pomocą funkcji `read.table`. Funkcja ta pozwala wczytać pliki, w których wartości oddzielone są białymi znakami (np. tabulatorami). Domyślnie jednak wczytuje pliki bez nagłówka (zakłada, że wszystkie wiersze to dane). My jednak ustawiliśmy argument `header` na `TRUE`, co sprawiło, że R potraktował pierwszy wiersz jako wiersz z nazwami kolumn.

```{r}
library(pander)
library(knitr)
data <- read.table("Tab15-1.dat", header = T)
kable(data)
```

Wczytane przez nas dane dotyczą edukacji w Stanach Zjednoczonych i zostały zebrane przez Gubner (1999). Jak widzimy, w każdym wierszu mamy informacje dotyczące jednego stanu.

Nas dzisiaj będą interesowały dane z kolumn `SATcombined`, `PTratio`, `PctSAT`, oraz `LogPctSAT`. 

W kolumnie `SATcombined` znajduje się średni wynik testu SAT w danym stanie. Test ten umożliwia zarekrutowanie się na studia wyższe. Nie jest to jednak jedyny taki test. Uczniowie mogą również zdawać test ACT. Generalnie rzecz biorąc test SAT jest popularniejszy w jednych stanach, w innych zaś większą popularnością cieszy się ACT. Różnice w popularności możemy ocenić przyglądając się wartościom z kolumny `PctSAT`. Znajduje się tam odsetek uczniów, którzy przystąpili do egzaminu SAT. W kolumnie `LogPctSAT` znajdziemy logarytm naturalny z tej wartości. Dodatkowo warto pamiętać, że SAT wymagany jest przez wiele prestiżowych uczelni. W kolumnie `PTratio` znajdziemy informacje o tym ilu uczniów przypada na jednego nauczyciela. W kolumnie `Expend` znajdują się wydatki na edukacje w danym stanie.

Zanim przystąpimy do analizy danych, obejrzyjmy rozkład naszych zmiennych. Dla każdej z nich narysujemy histogram, wykres kwantyl-kwantyl (który posłuży nam do oceny normalności rozkładu) oraz wykres rozrzutu, gdzie zmienną na osi Y będzie `SATcombined` (oprócz samego `SATcombined`).

Żeby to zrobić, posłużymy się odpowiednimi funkcjami z pakietu `ggplot2`. Jeżeli nigdy nie korzystali Państwo z tego pakietu, to zachęcam do przerobienia odpowiedniego notebooka ze strony kursu (nagrałem też do nich filmiki, które dostępne są na YT). Dla tych z Państwa, którzy korzystali z tego pakietu, kilka ciekawych twistów:

- Użyliśmy pakietu `gridExtra` i zawartej w nim funkcji `grid.arrange`, żeby narysować kilka wykresów stworzonych za pomocą `ggplot2` na jednym panelu. Jesli rysowalibyśmy te wykresy używając standardowych funkcji R takich jak `hist`, to moglibyśmy po prostu użyć `par(mfrow=c(4,3))`. W przypadku `ggplot2` technika ta jednak nie działa.
- Kazdy z wykresów jest obiektem, który można przypisać do zmiennej, co zrobiliśmy.
- Wszystkie te obiekty przekazaliśmy do funkcji `grid.arrange`, która narysowała te wykresy.

To, w jaki sposób rysuje się zawarte poniżej typy wykresów w `ggplot2`, omówione jest w notatniku poświęconym temu pakietowi, więc odsyłam do niego (oraz do dokumentacji `ggplot2`!), jeżeli ktoś chciałby dowiedzieć się więcej.

```{r}
library(ggplot2)
library(gridExtra)
expend_qq <- ggplot(data, aes(sample=Expend))+stat_qq() + stat_qq_line()
expend_hist <- ggplot(data, aes(x=Expend))+ geom_histogram(bins = 8)
expend_scatter <- ggplot(data, aes(x=Expend, y=SATcombined)) +geom_point() + geom_smooth(method=lm)


ptratio_qq <- ggplot(data, aes(sample=PTratio))+stat_qq() + stat_qq_line()
ptratio_hist <- ggplot(data, aes(x=PTratio))+ geom_histogram(bins = 8)
ptratio_scatter <- ggplot(data, aes(x=PTratio, y=SATcombined)) +geom_point() + geom_smooth(method=lm)

pctsat_qq <- ggplot(data, aes(sample=PctSAT))+stat_qq() + stat_qq_line()
pctsat_hist <- ggplot(data, aes(x=PctSAT))+ geom_histogram(bins = 8)
pctsat_scatter <- ggplot(data, aes(x=PctSAT, y=SATcombined)) +geom_point() + geom_smooth(method=lm)

satcombined_qq <- ggplot(data, aes(sample=SATcombined))+stat_qq() + stat_qq_line()
satcombined_hist <- ggplot(data, aes(x=SATcombined))+ geom_histogram(bins = 8)
logpctsat_scatter <- ggplot(data, aes(x=LogPctSAT, y=SATcombined)) +geom_point() + geom_smooth(method=lm)

grid.arrange(expend_hist, expend_qq, expend_scatter, # dajemy tutaj wszystkie wykresy
             ptratio_hist, ptratio_qq, ptratio_scatter,
             pctsat_hist, pctsat_qq, pctsat_scatter,
             satcombined_hist, satcombined_qq, logpctsat_scatter,
             ncol = 3 # liczba kolumn, w które chcemy je ułożyć
             )
```

Rzeczą, która powinna zwrócić naszą uwagę jest to, że relacja między `Expend` i `SATCombined` jest negatywna, to znaczy im więcej pieniędzy wydajemy na edukacje, tym gorsze wyniki uzyskują uczniowie (sic!). Jest to dość kontrintuicyjne, ponieważ spodziewalibysmy się zupełnie odwrotnego efektu -- wydatki na edukacje powinny przecież polepszać jej jakość.

Spróbumy teraz obliczyć współczynnik korelacji, między zmiennymi, które nas interesują. Jak Państwo wiedzą, do obliczenia współczynnika korelacji oraz przeprowadzania odpowiedniego testu statystyznego ($H_0: \rho = 0$) używamy funkcji `cor.test`. Dla przykładu możemy obliczyć korelację między zmienną `SATcombined` oraz `Expend`:

```{r}
pander(cor.test(~ SATcombined + Expend, data = data))
```

Widzimy, że uzyskaliśmy negatywny współczynnik korelacji --- niezbyt wysoki, ale statystycznie istotny! Jest to dość zaskakujące. 

Efekt ten, przynajmniej według podręcznika, z którego pochodzi ten przykład, związany jest z faktem, że w tych stanach, w których mały odsetek uczniów przystępuje do testu SAT, przystępują do niego zazwyczaj najlepsi i to zaburza wynik. W dalszej części notatnika pokażemy w jaki sposób możemy za pomocą wielokrotnej regresji to pokazać!

Krótka uwaga na marginesie. Gdybyśmy mieli policzyć współczynnik korelacji dla wszystkich kombinacji zmiennych, musielibyśmy napisać bardzo dużo kodu! Możemy sobie z tym jednak poradzić. W pakiecie `rstatix` znajduje się sporo ciekawych funkcji do ,,masowego'' przeprowadzania testów statystycznych. My skorzystamy z funkcji `cor_test` (która jest odpowiednikiem `cor.test` z biblioteki standardowej R).

```{r, message=F}
library(rstatix)
# Do funkcji cor_test przekazujemy ramkę danych oraz nazwy kolumn.
results <- cor_test(data = data, Expend, PTratio, Salary, PctSAT, SATcombined, LogPctSAT)
# Dane wyjściowe `cor_test` to ramka danych więc możemy ją sformatować za pomocą `kable`
kable(results[, c("var1", "var2", "cor", "p")])
```

Przejdźmy teraz do głównego tematu dzisiejszych zajęć, to znaczy do wielorakiej regresji liniowej. Jak już wspomnieliśmy, chcielibyśmy zobaczyć związek między wydatkami na szkolnictwo i testem SAT kontrolując przy tym odsetek przystępujących do SAT uczniów, to znaczy zmienną `LogPctSAT` (logarytm z odsetka przystępujących osób, wybraliśmy go bo ma troszkę lepszy rozkład i zalezność jest bardziej liniowa).

Stworzymy dwa modele w których zmienną objaśnianą jest `SATcombined`. W pierwszym z nich (`fit1`) uwzględniać będziemy tylko `Expend`.

```{r}
fit1 <- lm(SATcombined ~ Expend, data = data)
pander(summary(model1))
```

Teraz stworzymy drugi model, w którym to uwzględnimy jednak dodatkową zmienną - `LogPctSAT`.

```{r}
fit2 <- lm(SATcombined ~ Expend + LogPctSAT, data = data)
pander(summary(fit2))
```

Widzimy, że zmiana jest gigantyczna! W poprzednim wypadku współczynnik regresji był dla `Expand` ujemny ($-20.89$) ale teraz, kiedy kontrolujemy zmienną `LogPctSAT` jest delikatnie dodatni ($11.13$). To jest dokładnie to, czego oczekiwaliśm! Relacja między `LogPctSAT` a `SATcombined` jest negatywna - im mniej osób przystępuje w danym stanie do egzaminu, tym wyższy wynik średnio uzyskują. To również jest zgodne z naszymi oczekiwaniami. Warto zwrócić uwagę, że dramatycznie zwiększyło się również dopasowanie naszego modelu do danych - współczynnik $R^2$ wzrósł z $0.145$ do $0.886$.

### Interpretacja wielokrotnej regresji

Spróbujmy prześledzi, co tak naprawdę robi wielokrotna regresja liniowa po to, aby zrozumieć jak interpretować jej wyniki. 

Stwórzmy model, w którym zmienną objaśnianą będzie `SATcombined` a predyktorem `LogPctSAT`.

```{r}
fit_sat <- lm(SATcombined ~LogPctSAT, data = data)
```

Obliczmy jakie wartości dla `SATcombined` przewiduje nasz model. Posłużymy się tutaj funkcją `predict`. Działa ona tak, jakbyśmy po prostu wzięli wyraz wolny regresji oraz współczynniki kierunkowe a następnie dla każdej obserwacji obliczali przewidywaną wartość.

```{r}
predictsat <- predict(fit_sat)
```

Następnie obliczmy różnicę między prawdziwymi (rzeczywiście zaobserwowanymi) wartościami a wartościami przewidywanymi przez regresje (czyli residua/reszty!). Przypiszmy je do zmiennej `residsat`.

```{r}
residsat <- data$SATcombined - predictsat
```

W tym momencie musimy dokonać waznej obserwacji. Wartości naszej nowej zmiennej (`residstat`) są niezależne od wartości `LogPctSAT` (ponieważ w pewnym sensie usunęliśmy z nich wpływ `LogPctSAT`). Dokładnie to samo możemy zrobić z naszą drugą zmienną, czyli `Expend`.

```{r}
fit_expend <- lm(Expend ~ LogPctSAT, data = data)
predictexpend <- predict(fit_expend)
residexpend <- data$Expend - predictexpend
```

Mamy teraz dwie zmienne, które są niezależne od `LogPctSat`. W takim razie, możemy stworzyć model, w którym przewidujemy wartości `residsat` na podstawie wartości `residexpend`!

```{r}
fit_resid <- lm(residsat ~ residexpend)
pander(summary(fit_resid))
```

Jak widać współczynnik regresji który uzyskaliśmy dla tego modelu z jedną zmienną jest dokładnie taki sam jak wóœczas, gdy mieliśmy dwie zmienne. Dlaczego? Dlatego, że ,,manualnie'' zrobiliśmy to, co regresja robi za nas -- obliczyliśmy współczynnik regresji kontrolując efekty, które można przypisać trzeciej zmiennej (bo ,,usunęliśmy'' te efekty zarówno z predyktora jak i ze zmiennej objaśnianej).

O wielokrotnej regresji liniowej możemy też myśleć w inny sposób. Wiemy, że w przypadku naszego modelu z dwoma predyktorami równanie regresji ma postać

$$
\hat{Y} = b_1X_1 + b_2X_2 + b_0
$$
Czyli umując to w kategoriach naszych zmiennych:

$$
\widehat{SAT} = b_1{Expend} + b_2{LogPctSAT} + b_0
$$

W naszym wypadku znamy $b_1$ oraz $b_2$ (bo są to nasze współczynniki regresji) więc równanie to możemy przedstawić również tak:

$$
\widehat{SAT} = 11.129 \times {Expend} - 78.203 \times {LogPctSAT} + 1147.100
$$

Możemy więc korzystając z tego wzoru obliczyć $\widehat{SAT}$. 

```{r}
predictions <- predict(fit2)
```

Obliczmy więc korelację między $\widehat{SAT}$ (naszymi predykcjami w zmiennej `predictions`) oraz $SAT$ (czyli wartościami zaobserwowanymi). Tak obliczony współczynnik korelacji podnieśmy do kwadratu.

```{r}
cor.test(data$SATcombined, predictions)$estimate ** 2
```

Przypomnijmy sobie jak wyglądał nasz model dla dwóch zmiennych, który stworzyliśmy wczesniej.

```{r}
pander(summary(lm(SATcombined ~ Expend + LogPctSAT, data = data)))
```

Jak widzimy uzyskaliśmy dokładnie tę samą wartość co w przypadku $R^2$ wcześniej! Nakierowuje nas to na pewien sposób myślenia o regresji, zgodnie z którym zasadniczo sprowadza się ona do utworzenia nowej zmiennej ($\widehat{SAT}$), będącej najlepszą liniową kombinacją predyktorów w naszym modelu (`Expend` oraz `LogPctSAT`). Przez ,,najlepszą liniową kombinację'' rozumiemy to, że nie ma żadnych innych wag, które moglibyśmy przypisać naszym predyktorom (współczynników regresji), które lepiej przewidywałyby naszą zmienną objaśnianą.

O współczynniku $R$ dla wielokrotnej regresji liniowej możemy myśleć po prostu jako o współczynniku korelacji $r$ Pearsona tyle że dla zmiennej stanowiącej najlepszą liniową kombinację predyktorów.