---
title: "Regresja liniowa - ciąg dalszy"
output:
  pdf_document: default
  html_document: default
---

W tej części kursu zajmiemy się kilkoma bardziej zaawansowanymi aspektami prostej regresji liniowej z jednym predyktorem.

**Błąd standardowy estymacji** - odchylenie standardowe reszt. Wielkość ta informuje o przeciętnej wielkości odchyleń empirycznych wartości zmiennej zależnej od wartości wyliczonych z modelu (teoretycznych). Parametr ten jest ważny w analizie regresji, gdyż stanowi miarę rozproszenia elementów populacji wokół linii regresji. Na jego podstawie możemy ocenić stopień "dopasowania" modelu do danych empirycznych.

Dokładność oszacowania regresji mozna ocenić za pomocą współczynnika $\phi^2$. Punktem wyjścia jest wariancja Y. Możemy pokazać, że:

$$
S^2(Y) = S^2(\hat{Y}) + S^2(U)
$$

Równość ta (*równość wariancyjna*) głosi, że całkowity obszar zmienności zmiennej zaleznej jest sumą zmienności wyjaśnianej regresją i zmienności resztowej (niewyjaśnianej przez regresję).

**Współczynnikiem zbieżności** $\phi^2$ lub **indeterminacji** nazywamy stosunek tej części zmienności badanego zjawiska, która nie jest wyjaśniana przez zmiany zmiennych objaśniających w funkcji regresji do całkowitej zmienności zmiennej objaśnianej.

$$
\phi^2 = \frac{S^2(U)}{S^2(Y)}
$$

**Współczynnik determinacji** $R^2$ mierzy, jaka część ogólnej zmienności zmiennej zależnej jest wyjasniana przez regresję liniową. Obliczamy go według wzoru:

$$
R^2 = \frac{S^2(\hat{Y})}{S^2(Y)} = 1 - \phi^2
$$

**Interpretacja współczynnika $R^2$**

* Im większe $R^2$ tym lepiej. Musimy jednak pamiętac, że dołączenie nowej zmiennej do modelu zawsze spowoduje zwiększenie jego wartości. Dlatego w przypadku kilku predyktorów lepiej używać poprawionego $R^2$, który uwzględnia, że $R^2$ obliczony jest z próby i jest "zbyt dobrze" dopasowany.

$$R^{2}_{popr} = R^2 - \frac{1}{n-2}(1-R^2) $$

**Założenia regresji**

* **homogeniczność wariancji w grupach** - założenie że wariancja Y dla każdej wartości X jest stała w populacji
* **normalność w grupach** - założenie, że wartości Y dla każdej wartości X mają rozkład normalny
* (takie meta-założenie) **relacja między X a Y jest relacją, którą oddaje funkcja liniowa**

**Narzędzia diagnostyczne:**

* **wykres przedstawiajcy na jednej osi wartości dopasowane przez model, a na drugiej residua (reszty) lub standaryzowane residua**
   * dla modelu adekwatnego średnia wartość resuduum nie powinna zależeć od wartości dopasowania (powinniśmy w wyniku dostać pas punktów losowo rozmieszczonych wokół prostej $y = 0$
   * jeżeli można dostrzec jakiś wzór w ułożeniu punktów, może to oznaczać, że relacja między zmiennymi nie jest liniowa
   * w szczególności niepokoją nas regularne krzywe
   
* **wykresy kwantylowe dla standaryzowanych residuów**
    * powinny wskazywać na normalność
    * nie przejmujemy się za bardzo niewielkimi odchyleniami przy wartościach skrajnych
    * przejmujemy się, jeżeli wykres kwantyl-kwantyl reszt wskazuje na wyraźną skośność tzn. nienormalność

* **wykres, na którym dla każdej wartości zmiennej objaśniającej wyznaczono pierwiastek z wartości bezwzględnej jego residuum standaryzowanego**
    * nie powinniśmy zaobserwować żadnego trendu, jesli takowy występuje oznacza to, że wariancja błędu nie jest stała
    * powinniśmy zaobserwować równomiernie rozłożone punkty i poziomą linię przebiegającą przez nie
    
* **wykres dźwigni**
    * nie ma znaczenia wzór - interesują nas obserwacje odstające
    * przerywane linie reprezentują dystans Cooka - nie chcemy, żeby nasze obserwacje odstające się tam znalazły, bo to by znaczyło, że miały duzy wpływ na nasza regresję
    * wyniki regresji mogą znacznie się zmienić, jesli wykluczymy te obserwacje - czasami warto to zrobć

* **wykres pudełkowy/skrzypcowy Y warunkowanej na X** 
    * mozna użyc, jeśli X ma tylko kilka poziomów 
    * pudełka powinny wyglądać w miarę podobnie
    
Więcej o wykresach diagnostycznych można przeczytać tutaj w języku angielskim i w książkach też:

http://data.library.virginia.edu/diagnostic-plots/

## Przykład

Większość z nas potrafi całkiem nieźle sobie poradzić z egzaminem, jeżeli wcześniej uczylismy się materiału. Jak dobrze jednak poradzilibyśmy sobie, gdybyśmy przystąpili do egzaminu nawet nie patrząc wcześniej na materiał?

Katz, Lautenschlager, Blackburn i Harris (1990) postanowili odpowiedzieć na to pytanie i poprosili studentów o zapoznanie się z krótkim fragmentem tekstu a następnie o odpowiedzenie na dotyczące go pytania wielokrotnego wyboru. W drugiej grupie studentów autorzy również poprosili o odpowiedź na te pytanie, uczestnikom nie dali jednak wcześniej przeczytania rzeczonego fragmentu. Pytania na teście były bardzo podobne do tych, z którymi stykają się uczniowie w USA na teście SAT. Wyniki uzyskane przez badaczy doprowadziły ich do wniosku, że uczniowie, którzy dobrze poradzili sobie na SAT poradzą sobie równie dobrze na tym teście, ponieważ oba angażują podstawowe umiejętności rozwiązywania testów, takie jak eliminacja nieprawdopodobnych odpowiedzi itp.

Nas będą interesowały odpowiedzi tylko tych studentów, którzy nie czytali fragmentu tekstu. Zmienną objaśnianą będzie wynik na SAT, zmienną objaśniającą zaś - wynik na teście, który badacze dali studentom.

Na początek przyjrzyjmy się naszym danym - wczytajmy je do R, zobaczmy pierwsze 6 wierszy oraz obliczymy podstawowe statystyki deskryptywne.

```{r}
df <- read.table('Tab9-6.dat', header = TRUE)
head(df)
summary(df)
```

Nastepnie dopasujmy model liniowy do naszych danych:

```{r}
fit <- lm(SATV ~ Score, data= df)
summary(fit)
```

Zobaczmy, jak wygląda wykres punktowy (wykres rozrzutu) dla zmiennej objaśniającej i objaśnianej oraz dorysujmy linię regresji.

```{r}
plot(SATV ~ Score, data = df)
abline(fit)
```

Teraz użyć możemy naszych narzędzi diagnostycznych, które omówiliśmy wcześniej i zobaczyć, czy założenia regresji liniowej nie zostały poważnie naruszone.

```{r}
par(mfrow = c(2,2))
plot(fit)
```

Co możemy stwierdzić na podstawie stworzonych przez nas wykresów diagnostycznych?

## Predykcja

Modelu liniowego możemy użyć do wyprowadzania z danych przewidywań. Na początek spróbujmy obliczyć przewidywane wartości zmiennej objaśnianej ($\hat{Y}$) dla każdej naszych wartości zmiennej objaśniającej ($X$). R słusznie przypomni nam, że wartości te dotyczą przyszłych wartości.

```{r}
predict(fit, interval = "predict")
```

Za pomocą tej samej funkcji możemy również dowiedzieć się, jakie przewidywane wartości daje nasz model dowolnej innej wartości predyktora.

```{r}
predict(fit, data.frame(Score = 60), interval = "predict")
```

Rozważając predykcje w regresji liniowej, możemy skonstruować dwa przedziały ufności. Jeden z nich to przedział ufności predykcji dla poszczególnych wartości, drugi dla średniej. **UWAGA** - nie mylić z przedziałem ufności dla linii regresji.

O różnicy między tymi dwoma przedziałami więcej można przeczytać np. tu:

http://www.ma.utexas.edu/users/mks/statmistakes/CIvsPI.html.

Poniżej znajduje się wykres ilustrujący różnicę, między tymi dwoma pojęciami. Szerszy przedział to przedział dla predykcji, węższy - dla średniej.

```{r}
points <- seq(30, 70, length.out = 1000)
predictions <- predict(fit, data.frame(Score = points), interval = "predict")
plot(points, predictions[,1], type = 'l')
points(points, predictions[,2], type = 'l', lty = 2)
points(points, predictions[,3], type = 'l', lty = 2)
predictions = predict(fit, data.frame(Score = points), interval = "confidence")
points(points, predictions[,2], type = 'l', lty = 3)
points(points, predictions[,3], type = 'l', lty = 3)
```

To samo możemy zrobić za pomocą funkcji z pakietu `rms`. Bardzo polecam zapoznanie się w wolnym czasie z możliwościami tego pakietu - dużo różnych rodzajów regresji i przydatne funkcje.

```{r message = F}
library(rms)
dd <- datadist(df)
options(datadist = "dd")

model_rms <- ols(SATV ~ Score , data = df)
model_rms
```

```{r}
plot(Predict(model_rms))
plot(Predict(model_rms, conf.type = 'individual'))
ggplot(Predict(model_rms))
ggplot(Predict(model_rms, conf.type = 'individual'))
```