---
title: "Regresja liniowa - $R^2$ i diagnostyka"
author: "Bartosz Maćkiewicz"
date: "20 02 2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

## Narzędzia diagnostyczne

Spróbujmy użyć poznanych narzędzi diagnostycznych do danych, o których na pewno wiemy że nie spełniają załozeń regresji liniowej. Skąd wziąć takie dane? Najlepiej samemu je wygenerować. Zmienna X będzie pochodzić z rozkładu jednostajnego ciągłego ($a = 0$, $b = 10$). Nasze Y stworzymy w następujący sposób - dla każdego X podniesiemy go do potęgi drugiej a następnie dodamy do niego liczbę z rozkładu jednostajnego ciągłęgo z zakresu $a = 0$, $b = 5x$. Widzimy, że w tej sytuacji pogwałcone mamy (chyba) wszystkie założenia regresji liniowej, o których mówiliśmy:

- relacja nie jest liniowa (bo oddaje ją funkcja kwadratowa)
- rozkład Y ze względu na X nie jest normalny (bo jest z rozkładu jednostajnego ciągłego)
- wariancja w grupach nie jest homogeniczna (bo zależy od wartości X!)

```{r}
x = runif(100, 0,10)
y = x ** 2 + runif(100, 0, 5*x)
```

Nikt nam oczywiście nie broni dopasować linii regresji do takich danych. Zróbmy to więc i zobaczmy jak wygląda nasz model.

```{r}
m = lm(y ~ x)
summary(m)
plot(x, y)
abline(m)
```

Okazało się, że zależnośc jest statystycznie istotna! Czy powinno nas to uprawnić do konkluzji, że między X i Y zachodzi liniowa zależność? Niekoniecznie. Spójrzmy co powiedza na ten temat wykresy diagnostyczne:

```{r, fig.height=10, fig.width=10}
par(mfrow = c(2,2))
plot(m)
```

Wszystkie wykresy diagnostyczne wyglądają nie tak jak jak powinny! Widać wyraźną tendencję na wykresie "Residuals vs Fitted" (obserwacje układają się w literę U), wykres "Normal Q-Q" pokazuje, że standaryzowane reszty nie mają rozkładu normalnego. Na pozostałych dwóch wykresach również widać wyraźną tendencję, której nie powinno tam być. Takie wykresy powinny nas nakierować na to, że relacja między zmiennymi nie jest liniowa. 

Moglibyśmy sobie (do pewnego stopnia) poradzić z tym problemem po prostu wyciągając pierwiastek z y:

```{r}
m = lm(sqrt(y) ~ x)
summary(m)
plot(x, sqrt(y))
abline(m)
```

Wykresy diagnostyczne pokazują, że decyzja przyniosła oczekiwany przez nas skutek, chociaż wciąż nie jest idealnie:

```{r, fig.height=10, fig.width=10}
par(mfrow = c(2,2))
plot(m)
```

## $R^2$

Poniżej krótka ilustracja o co chodzi w $R^2$. Załóżmy, że mamy dane z omawianego eksperymentu dotyczącego *guilty pleasure*:

```{r}
data = read.csv("Study3_Results_AfterExclusion.csv")
data$IDEAL = apply(data[,c("IDEAL1", "IDEAL2")], 1, mean)
model = lm(FEELBAD ~ IDEAL, data = data)
plot(jitter(FEELBAD) ~ jitter(IDEAL), data = data)
abline(model)
summary(model)
```

Widzimy, że $R^2$ wynosi w naszym przypadków $0.23$ co interpretować możemy tak, że IDEAL wyjaśnia 23% zmienności albo wariancji FEELBAD. Co to jednak dokładnie znaczy? Obliczmy sumę kwadratów reszt dla FEELBAD bez uwzględniania naszej zmiennej IDEAL:

```{r}
SS_Y = sum((data$FEELBAD - mean(data$FEELBAD))**2)
SS_Y
```

To jest dokładnie tak, jakby Państwo obliczali wariancję dla FEELBAD (tylko bez dzielenia).

Teraz zastanówmy się jaka jest tak obliczana wariancja jeśli:
- weźmiemy pod uwagę tylko wartości przewidywane z naszego równania regresji (SS_Y_P jak *predicted*)
- weźmiemy pod uwagę tylko reszty z regresji (SS_res jak *residuum*)

```{r}
predicted = predict(model, data) # w predicted są wartości Y tak jak je przewiduje nasze równanie regresji
SS_Y_P = sum((predicted - mean(data$FEELBAD))**2)
SS_Y_P

SS_res = sum((data$FEELBAD - predicted)**2)
SS_res
```

Łatwo zauwazyć, że suma tych dwóch sum kwadratów reszt jest dokładnie taka sama jak nasza wyjściowa wariancja:

```{r}
SS_Y_P + SS_res
```

Innymi słowy rozłożyliśmy naszą wariancje na wariancję, którą można przypisać zmiennej IDEAL oraz pozostałą wariancję. $R^2$ jest więc stosunkiem wariancji wyjaśnianej przez predyktor w regresji ą i całej wariancji. Całkowita wariancja zmiennej ($SS_Y$) zawsze pozostaje taka sama, ale im większe jest nasze $SS_{\hat{Y}}$ i tym samym mniejsze nasze $SS_{res}$ tym nasz model jest lepiej dopasowany do danych!

```{r}
r2 = SS_Y_P / SS_Y
r2
```

Otrzymaliśmy (tak jak się spodziewaliśmy!) tę samą wartość, którą obliczył nam R.

## Adjusted r ($r_{adj}$)

To jest nieobciążony estymator korelacji w populacji (którego z jakiegoś powodu się nie używa, ale w moich tajnych książkach go znalazłem):

$$
r_{adj} = \sqrt{1 - \frac{(1-r^2)(N-1)}{N-2}}
$$

A to jest $R^2_{adj}$ ($p$ oznacza liczbę predyktorów):
$$
R^2_{adj} = 1 - (1- R^2)\frac{n-1}{n-p-1}
$$
dla jednego predyktora jest to:

$$
R^2_{adj} = 1 - (1- R^2)\frac{n-1}{n-2}
$$
Kiedy wie się, jak wygląda nieobciążony estymator $\rho$  (czyli korelacji w populacji), to nie powinno dziwić :)