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

W tym odcinku kontynuować będziemy nasze przygody z wielokrotną regresją liniową. Tym razem wejdziemy troszkę głębiej w diagnostykę, założenia regresji oraz interpretację wyników. Cały czas korzystać będziemy z naszych danych dotyczących edukacji, które ściągnąć można z platformy. Na początek wczytajmy więc nasze dane i dla przypomnienia wyświetlmy sobie pięć pierwszych wierszy. Osoby, które nie pamiętają już czego dane dotyczą proszone są o zajrzenie do podręcznika i zajrzenie do notatnika z poprzedniego tygodnia.

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

## Tolerancja i VIF

Współczynnik inflacji wariancji (*Variance Inflation Factor*) pozwala nam ocenić, czy predyktory są *współliniowe*. To, w jaki sposób ten współczynnik obliczyć nie jest w tej chwili istotne, potraktujmy go po prostu jako pewną wartość diagnostyczną. Jeżeli między predyktorami w naszym modelu zaobserwować możemy  duży stopień współliniowości (powiedzmy, żę $VIF > 4$), to wówczas może to skutkować niestabilnością naszego modelu. 
Tolerancja jest odwrotnością VIF i oblicza się ją $1/VIF$. Co do zasady dążymy do tego, aby tolerancja była jak największa a VIF jak najmniejszy. Mówi nam to, że nasz model charakteryzuje wysoka stabilność (niski błąd standardowy współczynników regresji). W przypadku $VIF > 10$ (silna współliniowość) powinniśmy dokonać korekty naszego modelu i usunąć ze zbioru predyktorów jedną zmienną.

VIF i Tolerancja pozwalają nam ocenić, które wybrane przez nas predyktory pokrywają się i w jakim stopniu. 

W R VIF i tolerancje obliczyć może np. za pomocą funkcji `vif` z pakietu `car`. Stwórzmy model, w którym objaśniamy `SATcombined` za pomocą `Expend`, `LogPctSAT` oraz `PTratio`.

```{r}
library(car)

fit <- lm(SATcombined ~ Expend + LogPctSAT + PTratio, data = data)
pander(summary(fit))
```

Następnie obliczymy współczynnik inflacji wariancji (VIF).

```{r}
pander(vif(fit))
```

Oraz Tolerancję (podzielimy po prostu 1 przez uzyskane przez nas współczynniki VIF).

```{r}
pander(1/vif(fit))
```

Jak widzimy VIF dla naszych predyktorów nie jest specjalnie wysoki, co uprawnia nas do stwierdzenie, że predyktory nie są ze sobą skorelowane w sposób zagrażający rzetelności naszego modelu.

## Standaryzowane współczynniki regresji i ,,ważność'' poszczególnych predyktorów

Jednym z istotniejszych pytań, jakie zadać sobie możemy przeprowadzając analizę regresji, jest pytanie o to, które z wybranych przez nas predyktorów są najważniejsze. Pytanie to oczywiście interpretować można na wiele sposobów. 

Błędem byłoby tutaj patrzeć po prostu na współczynnik regresji i uznać, że te predyktory, które mają wyższy współczynnik regresji, są ważniejsze. Jest to niewłaściwe rozumowanie, ponieważ współczynniki regresji zależą od tego w jakich jednostkach wyrażone są wartości naszych predyktorów. Możemy jednak sobie z tym  poradzić standaryzując zmienne, których używamy w charakterze predyktorów (przypominam: odejmujemy od wszystkich wartości średnią i dzielimy wszystkie wartości przez odchylenie standardowe dzięki czemu nasze zmienne mają średnią 0 i odchylenie standardowe 1). Dodatkowo wystardaryzować możemy zmienną objaśnianą. Tak wystandaryzowane zmienne posłużyć nam mogą do przeprowadzenia regresji.

```{r}
data$SATcombinedStd <- (data$SATcombined - mean(data$SATcombined))/sd(data$SATcombined)
data$ExpendStd <- (data$Expend - mean(data$Expend))/sd(data$Expend)
data$LogPctSATStd <- (data$LogPctSAT - mean(data$LogPctSAT))/sd(data$LogPctSAT)
data$PTratioStd <- (data$PTratio - mean(data$PTratio))/sd(data$PTratio)

fit_std <- lm(SATcombinedStd ~ ExpendStd + LogPctSATStd + PTratioStd, data = data)
pander(summary(fit_std))
```

Jak widzimy uzyskaliśmy trochę inne współczynniki regresji. Tym razem mają one postać standaryzowanych współczynników regresji. Warto zwrócić uwagę na to, że wyraz wolny regresji (*intercept*) wynosi w tym wypadku 0. Dlaczego? Dlatego, że wystandaryzowaliśmy zmienną objasnianą. Dodatkowo zmienia się teraz interpretacja wartości w kolumnie `Estimate` regresji - teraz, kiedy mamy tam standaryzowane współczynniki regresji, mówi nam ona, o ile odchyleń standardowych zmienia się nasza zmienna objaśniania wraz ze zmianą naszych predyktorów również wyrażonej w odchyleniach standardowych. Teraz, kiedy wszystko wyrazone jest już w odchyleniach standardowych, możemy sensownie porównywać predyktory między sobą.

Proszę zwrócić uwagę również na to, że standaryzacja zmiennych w regresji nie ma wpływu wartości statystyki testowej $t$ ani na wartości $p$, chociaż oczywiście ma wpływ na błąd standardowy (bo zmieniliśmy jednostki!).

W zdefiniowanym przez nas sensie najważniejszym predyktorem jest `LogPctSAT` ponieważ jego wzrost o jedno odchylenie standardowe skutkuje zmianą w `SATcombined` wynoszącą $-1.042$ odchylenia standardowego.

To samo co zrobiliśmy ,,na piechotę'' w R możemy również zrobić za pomocą funkcji `lm.beta` z pakietu `lm.beta`. Jako argument przekazujemy nasz modela a funkcja ta sama uzupełni go o standaryzowane współczynniki regresji.

```{r, warning=F}
library(lm.beta)
pander(summary(lm.beta(fit)))
```

## Błędy standardowe współczynników regresji i testy statystyczne dla współczynników regresji

Przyjrzyjmy się ponownie stworzonemu przez nas modelowi.

```{r}
pander(summary(lm.beta(fit)))
```

Widzimy, że w tabeli oprócz `Estimate` czyli współczynników regresji mamy również `Std.Error`, `t value` oraz `Pr(>|t|)`. Jak Państwo pamiętają `Pr(>|t|)` mówi nam, czy dany współczynnik regresji różni się w statystycznie istotny sposób od $0$. R oblicza go za nas, ale nie jest trudno przeprowadzić taki test samemu. `Std.Error` jest błędem standardowym współczynnika kierunkowego regresji. Proszę mysleć o tej wartości tak samo, jak myślą Państwo o błędzie standardowym średniej -- jest to (mówiąc trochę obrazowo) wartość odchylenia standardowego dla danej statystyki gdybyśmy brali taką samą próbę nieskończenie wiele razy i za każdym razem obliczali daną statystykę. Z konceptualnego punktu widzenia różnica jest taka, że w znanym Państwu przypadku było to odchylenie standardowe średnich z próby, tutaj jest to odchylenie standardowe współczynników kierunkowych regresji.

Statystykę testową $t$ obliczyć możemy według wzoru:

$$
t = \frac{b_j - b_{j}^*}{s_{b_j}}
$$

gdzie $b_j$ jest wartością naszego współczynnika kierunkowego regresji, $b_{j}^*$ jest wartością przy założeniu hipotezy zerowej a $s_{b_j}$ błędem standardowym. Dla $H0: b_{j}^*  = 0$ (czyli dla hipotezy zerowej głoszącej, że współczynnik kierunkowy regresji wynosi 0 czyli nie ma relacji między predyktorem a zmienną objaśnianą) wzór ten ma postać:

$$
t = \frac{b_j}{s_{b_j}}
$$

Rozkład statystyki testowej $t$ przy założeniu hipotezy zerowej ma w tym wypadku $N - p - 1$ stopni swobody gdzie $N$ to liczba obserwacji a $p$ to liczba predyktorów.

Sprawdźmy to! Wiemy, że dla zmiennej `Expend` współczynnik kierunkowy wynosi $11.67$ a błąd standardowy $3.533$. Obliczmy więc że statystykę testową:

```{r}
tstat = 11.67/3.533
tstat
```

Teraz możemy sprawdzić za pomocą funkcji `pt` jakie jest prawdopodobieństwo uzyskania takiej lub bardziej skrajnej wartości statystyki testowej: 

```{r}
pt(tstat, nrow(data) - 3 - 1, lower.tail = F) * 2 # mnożymy przez dwa bo dwustronny test
```

Jak widzimy nasz wynik zgadza się z tym, który obliczył dla nas R. Uzyskaliśmy wartość $p$ mniejszą niż $0.05$, co uprawnia nas (przy założonym poziomie istotności statystycznej $\alpha = 0.05$) do odrzucenia hipotezy zerowej.

## Założenia regresji liniowej

Jeżeli nasze predyktory są zmiennymi losowymi, to zakładamy, że łączny rozkład $Y, X_1, X_2, ..., X_p$ jest *wielowymiarowym rozkładem normalnym* (*multivariate normal distribution*).

Jeżeli poziomy naszych predyktorów są ustalone (czyli np. manipujemy nimi w ramach schematu eksperymentalnego), to wówczas musimy założyć, że rozkład warunkowy dla każdego poziomu  predyktorów jest w naszej zmiennej objaśnianej Y normalny.

W praktyce jednak nawet dość spore odchylenia od tego założenia nie wpływają negatywnie na poziom błędu I rodzaju dla przeprowadzanych przez nas testów statystycznych. Warto jednak o tym pamiętać!

Na wikipedii można znaleźć ładne obrazki pokazujące jak wyglądają takie rozkłady. Zachęcam do pogooglowania (sam nie podejmuje się narysowania ich w R ponieważ są trójwymiarowe...)

## Współczynnik $R^2$ i testy statystyczne

O $R^2$ mówiliśmy już sporo w poprzednich cześciach kursu, dlatego odsyłam do materiałów sprzed tygodnia i dwóch. Teraz jednak ponownie przyjrzyjmy się stworzonemu przez nas modelowi.

```{r}
pander(summary(lm.beta(fit)))
```

Widzimy, że R podaje nam dwie wartości statystyki $R^2$ - zwykłąd oraz ,,adjusted''. Wynika to z tego, że $R^2$ nie jest nieobciążonym estymatorem odpowiedniego parametru w populacji. Wersję nieobciążoną tego współczynnika oblicza się według wzoru:

$$
R^2_{adj} = 1 - \frac{(1-R^2)(N-1)}{N - p - 1}
$$

gdzie $N$ to liczba obserwacji a $p$ to liczba predyktorów. Z jakiegoś powodu jednak w artykułach naukowych wartość ta pojawia się dość rzadko (tutaj zgadzam się z Howellem).

Możemy także przeprowadzić test statystyczny dla uzyskanego przez nas współczynnika $R^2$. R robi to za nas i możemy to zobaczyć jeśli wyświetlimy sobie pełne `summary` naszego modelu (pakiet `pander` trochę ucina, niestety).

```{r}
summary(fit)
```

Naszą hipoteza zerowa głosi w tym przypadku, że $H_0: R^2 = 0$ (w tym wypadku $R^2$ odnosi się do tego współczynnika w populacji). Odpowiednią statystykę testową obliczyć możemy według wzoru:

$$
F = \frac{(N-p-1)R^2}{p(1-R^2)}
$$
I (przy założeniu $H_0$) ma ona rozkład $F$ o $p$ i $N-p-1$ stopniach swobody. Warto tutaj zauważyć, że statystyczna istotność testu zależy tu nie wyłącznie od uzyskanego w próbie $R^2$ oraz liczebności próby, ale także od liczby predyktorów! Im mamy ich więcej, tym trudniej uzyskać nam statystycznie istotny wynik.

Zobaczmy, czy uda nam się uzyskać taki sam wynik jak wyliczył nam R.

```{r}
fstat <- ((nrow(data) - 3 - 1) * 0.8865) / (3 * (1 - 0.8865))
fstat
```

I teraz sprawdźmy istotność statystyczną dla obliczonej statystyki testowej $F$ :

```{r}
pf(fstat, 3, nrow(data) - 3 - 1, lower.tail = F) * 2
```

Uzyskaliśmy wynik mniejszy niż 0.05, co daje nam podstawy do odrzucenia hipotezy zerowej.

## Wielkość próbki

Jak dużą powinniśmy mieć próbę, kiedy w naszym badaniu chcemy zastosować regresję liniową? Odpowiedż na to pytanie nie jest jednoznaczna, ale w literaturze proponuje się kilka ,,reguł kciuka'', które mogą ułatwić podjęcie decyzji o planowanej liczebności próby.

- co najmniej 10 obserwacji dla każdego predyktora ($N \geq 10p$)
- liczba obserwacji powinna przekraczać liczbę predyktorów o przynajmniej 50 ($N \geq p + 50$)

Możemy również zastosować inne podejście i rozważyć to od strony mocy statystycznej. 

(Howell odwołuje się w podręczniku do książki Cohena, Cohen Westa i Aikena, ale nie byłem w stanie znaleźć podanych przez Howella liczb i mam podejrzenie, że je zmyślił. Znalazłem jednak inne przykłady  w tej książce, które umiem zreplikować w R. Jeżeli ktoś wie skąd Howell wziął te liczby, to będzie zwolniony z jednej pracy domowej)

Jak dużej potrzebujemy próby, aby mieć 80% szans na wykryce prawdziwego $R^2$ wynoszącego 0.2 przy trzech predyktorach? Odpowiedź na to pytanie da nam analiza mocy statystycznej. W R możemy to zrobić za pomocą dostępnej w pakiecie `pwr` funkcji `pwr.f2.test`. Funkcja ta przyjmuje szereg argumentów, ale trzy z nich mogą być na pierwszy rzut oka mylące.

- `u` - liczba stopni swobody w liczniku, to jest liczba predyktorów, które planujemy
- `v`- liczba stopni swobody w liczniku, to jest nasze $N - p - 1$ (proszę spojrzeć na wzór na statystykę testową F)
- `f2` - to jest wielkość efektu, którą oblicza się według wzoru $\frac{R^2}{1-R^2}$ (znów, prosze spojrzeć na wzór)

Teraz, kiedy już wiemy co oznaczają poszczególne argumenty tej funkcji, możemy z niej skorzystać. Musimy dokładnie jeden argument tej funkcji pozostawić pustym, my chcemy się dowiedzieć czegoś o wielkości próby, więc pustym pozostawiamy `v`.

```{r}
library(pwr)
pwr.f2.test(u = 3, v = NULL, f2 = 0.2/ (1 -0.2), power = 0.8, sig.level = 0.05)
```

Otrzymaliśmy `v` wynoszące 44 (po zaokrągleniu do góry) co pozwala nam łatwo obliczyć niezbędną liczebność próby: $44 + 3 + 1 = 48$.

W jaki sposób liczba predyktorów w naszym modelu wpływa na moc testu? Rozważmy dwa identyczne przypadki. Interesuje nas wykryce z prawdopodobieństwiem $0.8$ ,,prawdziwego'' $R^2$ wynoszącego przynajmniej 0.09. Ile obserwacji potrzebujemy dla jednego predyktora?

```{r}
pwr.f2.test(u = 1, v = NULL, f2 = (0.3)**2/ (1 -(0.3**2)), power = 0.8, sig.level = 0.05)
```

$80 + 1 + 1$ czyli $82$. Czy liczebność wzrośnie, kiedy założymy, że będziemy mieć nie jeden ale pięć predyktorów?

```{r}
pwr.f2.test(u = 5, v = NULL, f2 = (0.3)**2/ (1 -(0.3**2)), power = 0.8, sig.level = 0.05)
```

Jak widzimy teraz potrzebujemy $130 + 5 + 1$ czyli $136$ obserwacji. Widzimy więc, że im więcej predyktorów, tym większej próby potrzebujemy, żeby zachować taką samą moc. W przypadku regresji najlepszą ,,regułą kciuka'' jest jednak ,,im więcej, tym lepiej'' (ale też bez przesady!).

## Cząstkowa i półcząstkowa korelacja

Korelacja cząstkowa to korelacja między dwiema zmiennymi ($X$ i $Y$) przy usuniętym wpływie jednej lub większej liczby zmiennych ($Z$). Jest to współczynnik korelacji Pearsona pomiędzy składnikami resztowymi (resztami) $X_{res}$ i $Y_{res}$ z modeli regresji liniowej prognozujących odpowiednio $X$ za pomocą $Z$ oraz $Y$ za pomocą $Z$. 

Proszę zwrócić uwagę, że robiliśmy już to (obliczaliśmy korelację między resztami) w poprzednim notatniku, kiedy zajmowaliśmy się podstawami wielokrotnej regresji, dlatego nie będę powtarzał tutaj kodu i odsyłam do poprzedniego notatnika. Korelacja cząstkowa jest o tyle przydatna, że pozwala nam np. stwierdzić, czy, ,,oczyszczając'' nasze obserwacje z wpływu jakiejś trzeciej zmiennej, zachodzi relacja między dwiema zmiennymi. 

Na początek przypomnijmy sobie jak wyglądał model z dwiema zmiennymi:

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

W R możemy łatwo obliczyć korelację cząstkową (i od razu przeprowadzić test!) za pomocą funkcji `pcor.test` z pakietu `ppcor` (w kolumnie `estimate` znajdziemy współczynnik korelacji): 

```{r, warning=F, message=F}
library(ppcor)
pcor.test(data$SATcombined, data$Expend, data$LogPctSAT)
```

Widzimy, że korelacja między `SATcombined` a `Expand` jest statystycznie istotna przy usuniętym wpływie `LogPctSAT`. Możemy to również obliczyć w drugą stronę:

```{r}
pcor.test(data$SATcombined, data$LogPctSAT, data$Expend)
```

Również uzyskaliśmy statystycznie istotny wynik, ale współczynnik korelacji jest znacznie wyższy. 

Innym typem korelacji jest korelacja semicząstkowa. Największą różnicą między korelacją cząstkową a semicząstkową jest to, że w przypadku tej pierwszej usuwamy wpływ trzeciej zmiennej zarówno z predyktora jak i zmiennej objaśnianej, podczas gry w przypadku tej drugiej jedynie z predyktora. 

Jeżeli więc obliczymy korelację semiczęstkową między `SATcombined` i `Expend` kontrolując `LogPctSAT`, to jest to korelacja między `SATcombined` a tą częścią `Expend`, która jest niezależna od `LogPctSAT`.

```{r}
spcor.test(data$SATcombined, data$Expend, data$LogPctSAT)
```

W jaki sposób interpretować wielkość tych współczynników? W przypadku cząstkowej korelacji, jeżeli podniesiemy do kwadratu współczynnik *r*, to uzyskamy:

```{r}
pcor.test(data$SATcombined, data$Expend, data$LogPctSAT)$estimate **2
```

Oznacza to, że 20% zmienności w wynikach SAT, której nie da się wyjaśnić za pomocą logarytmu z odsetka osób podchodzacych do testów w danym stanie (`LogPctSAT`), może być wyjaśniona przez tę część nakładów finansowych na edukację (`Expend`), której nie da się wyjaśnić za pomocą logarytmu z odsetka osó” podchodzacych do testów w danym stanie.

Warto zwrócić uwagę, że w kontekście modelu liniowego, test istotności statystycznej dla korelacji cząstkowej dają taki sam wynik, jak testy dla współczynników regresji. Jest to nieprzypadkowe, jeżeli przypomnimy sobie poprzedni notatnik w którym właśnie używając procedury obliczania korelacji cząstkowej tłumaczyliśmy sobie jak działa wielokrotna regresja!

## Porównywanie modeli

Jednym z najczęściej przydarzających się problemów związanych z analizą regresji jest to, który z alternatywnych modeli powinniśmy wybrać. Możemy na przykład mieć dwie teorie które mają odmienne predykcje dotyczące tego, które zmienne będą odpowiadać za zmienną objaśnianą. Możemy oczywiście po prostu porównywać $R^2$ obu modelu, ale wiemy, że co do zasady zwiększenie liczby predyktorów zwiększy również dopasowanie modelu.

Rozważmy ten problem w dwóch sytuacjach. Pierwsza z nich to sytuacja, w której mamy dwa modele, ale zbiór predyktorów jednego z nich zawiera się w zbiorze predyktorów drugie. Taką sytuację nazywać będziemy modelami zagnieżdżonymi. Druga możliwośc jest taka, że zbiór predyktorów jednego nie zawiera się w zbiorze predyktorów drugiego.

### Modele zagnieżdżone (hierarchiczne)

W tym wypadku sytuacja jest relatywnie prosta. Rozważmy dwa modele. W pierwszym z nich w zbiorze predyktorów uwzględnimy zmienne `Salary`, `PTratio`, `LogPctSAT` oraz `Expend`. W drugim zaś będziemy mieli tylko `LogPctSAT` oraz `Expend`. Jak widzimy zbiór predyktorów drugiego modelu zawiera się w zbiorze predyktorów pierwszego, więc mamy do czynienia z zagnieżdżonymi modelami.

Pierwszy z nich wygląda tak:

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

Drugi zaś tak:

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

Widzimy, że pierwszy z nich ma wyższy współczynnik $R^2$. Czy jednak różnica ta jest statystycznie istotna? Pominiemy tutaj (dośc wyjątkowo jak na dzisiejszy notebook) przeprowadzanie odpowiedniego testu na piechote i skorzystamy z funkcji `anova` pochodzacej z biblioteki standardowej. Jako argumenty przyjmuje ona jeden lub więcej obiektów z modelami.

```{r}
anova(fit1, fit2)
```

Widzimy, że otrzymana przez nas wartość $p$ nie przekroczyła progu statystycznej istotności. Oznacza to, że w tym wypadku dodanie dodatkowych predyktorów nie zwiększyło istotnie dopasowania naszego modelu.

### Kryterium informacyjne Akaikego (AIC, Aikake's Information Criterion)

Co jednak zrobić, kiedy nie mamy do czynienia z modelami zagnieżdżonymi? Tutaj z pomocą przychodzą inne techniki porównywania ze sobą modeli. Jednym z narzędzi, którego możemy użyć jest kryterium informacyjne Akaikego. Rozważmy sytuację, w której mamy dwa modele - różniące się tym, czy użyliśmy zmiennej `LogPctSAT` czy `PctSAT` (logarytm z odsetka czy sam odsetek).

Możemy obliczyć wartość AIC za pomocą funkcji AIC.

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

pander(AIC(model1, model2))
```

Im mniejsza wartość AIC tym lepiej. Widzimy, że w naszym wypadku wersja z `LogPctSAT` ma mniejszy AIC niż wersja z `PctSAT`, co daje nam powody sądzić, że to właśnie ten model jest lepiej dopasowany.

(Jeżeli ktoś odkryje dlaczego SPSS Howella zwraca inną wartość niż R, to również ma odpuszczoną pracę domową).