--- 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.