---
title: "Regresja liniowa - powtórka"
author: "Bartosz Maćkiewicz"
date: "20 02 2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

## Przykład: Badania dotyczące *guilty pleasure*

Goffin i Cova w swoim artykule (2019!) zmierzyli się z nurtującym wszystkich pytaniem estetycznego nihilizmu. Dotychczasowe badania empiryczne wskazują, że ludzie nie uważają sądów estetycznych za ,,poprawne'' lub ,,niepoprawne'', ,,dobre'' lub ,,złe''. Nasze oceny dzieł sztuki nie odzwierciedlają żadnych ,,obiektywnych'' wartości estetycznych w nich zawartych.

Istnieje jednak fenomen *guilty pleasure* polegający na tym, że ktoś czerpie przyjemność z dzieł sztuki (szeroko rozumianych!) i jednocześnie czuje się z tego powodu źle. Dlaczego? Jedna hipoteza (*estetyczna*) głosi, że jest tak, ponieważ osoby te zdają sobie sprawę, że dziełom, z których czerpią przyjemność, brakuje pewnych obiektywnych wartości estetycznych i dlatego ich zachowanie jest nieodpowiednie. Druga (*nieestetyczna*) głosi zaś, że chodzi nie o transgresję norm estetycznych, ale osobistych ideałów lub norm społecznych, brak wartości intelektualnej dzieła czy odstawanie od poglądów innych ludzi. Na przykład przyznanie się, że słucha się Bonusa RPK przed innymi może doprowadzić do wysnucia wielu niekorzystnych wniosków na nasz temat.

Badacze przeprowadzili serię eksperymentów, w których próbowali rozstrzygnąć tę kwestię. W jednym z nich chcieli zobaczyć, które czynniki najlepiej wyjaśniają ,,czucie się źle'' związane z partycypacją w *guilty pleasures*.

Badani otrzymali od Goffina i Cova'y definicję *guilty pleasure* oraz poprosili badanych o przypomnienie sobie i opisanie jakiegoś swojego *guilty pleaser* oraz opisanie własnymi słowami, dlaczego czuli się z tym źle. Następnie badani na skali 1 (TOTALLY DISAGREE) do 7 (TOTALLY AGREE) ocenili 18 stwierdzeń dotyczących swojego obiektu *guilty pleasure*. Następujące dwa stwierdzenia miały charakter kontrolny i miały potwierdzić, czy faktycznie mamy do czynienia z tym zjawiskiem:

> (C1) "I really enjoy this work."

> (C2) "I feel bad about enjoying this work."

Część pytań dotyczyła różnych aspektów tego doświadczenia. Niektóre z nich dotyczyły jego normatywnego wymiaru jak na przykład:

> (N2) "I feel like my enjoyment of this work is wrong and musguided."

Inne dotyczyły identyfikacji z tym negatywnym doświadczeniem:

> (I1) "My enjoyment is in conflict with my own opinions and thoughts about the work."

Pozostałe pytania miały na celu przetestować różne wyjaśnienia zjawiska *guilty pleasure*. Wyjaśnienie estetyczne testowane było między innymi przez pytanie:

> (Art2) "If I were to write an art review. I would say the object of my enjoyment is a poor work of art."

Jedno z pytań testujące wyjaśnienie, zgodnie z którym *guilty pleasure* wynika z braku intersubiektywnej zgody:

> (Int2) "My enjoyment of this work does not match other people's experience."

Wyjaśnienie, które *guilty pleasure* interpretuje jako zaprzeczenie swoim osobistym ideałom tesotwane było przez pytania, wśród których znalazło się:

> (Personal1) "People I admire would probably not enjoy this kind of work."

Jeśli chodzi o wyjaśnienie biorące pod uwagę normy społeczne, to jedno z pytań brzmiało:

> (Social2) "I'm afraid of what people think of people who enjoy this kind of artwork."

Ostatnim testowanym wyjaśnieniem było to, zgodnie z którym dziełom tym brakuje wartości intelektualnych (czyli niekoniecznie estetycznych):

> (Intellectual1) "This work is not very complex or intellectual."

Badacze udostępnili cały komplet danych, więc skorzystamy z okazji i spróbujemy własnoręcznie przeanalizować te dane. Na początek wczytamy je do R:

```{r}
data <- read.csv("Study3_Results_AfterExclusion.csv")
```

### Regresja liniowa

Regresje przeprowadzamy w R za pomocą funkcji `lm`. Funkcja ta zwraca obiekt, którego metoda `print` wyświetla wyraz wolny regresji, współczynnik kierunkowy oraz wywołanie funkcji. Jeżeli chcemy się dowiedzieć czegoś więcej, musimy obiekt ten przekazać funkcji `summary`.

Na początek sprawdźmy, jak wygląda relacja między odpowiedziami na pytania testujące wyjaśnienie estetyczne oraz poziomem zgody na stwierdzenie o tym, że badani czuli się źle ze swoim *guilty pleasure*.


```{r}
data$ART <- apply(data[,c("ART1", "ART2", "ART3")], 1, mean)
data$IDEAL <- apply(data[,c("IDEAL1", "IDEAL2")], 1, mean)
data$PEOPLE <- apply(data[,c("PEOPLE1", "PEOPLE2")], 1, mean)
data$SOCIAL <- apply(data[,c("SOCIAL1", "SOCIAL2")], 1, mean)
data$INTELLECTUAL <- apply(data[,c("INTELLECTUAL1", "INTELLECTUAL2")], 1, mean)
lm(data$FEELBAD ~ data$ART)
```

**Wywołanie funkcji `summary` na obiekcie zwracanym przez `lm`**

```{r}
# Funkcja `summary` pozwala nam się dowiedzieć więcej o modelu zwróconym przez funkcję `lm`,
# w szczególności zaś poznać różne jego parametry
summary(lm(data$FEELBAD ~ data$ART))
```

### Dodawanie linii regresji

Moglibyśmy ręcznie wydobyć wartości niezbędne do dodania linii regresji, ale możemy skorzystać z faktu, że nasz model możemy bezpośrednio przezkazać funkcji `abline`. Należy jednak pamiętać o tym, żeby tworząc wykres punktowy używać składni formuły (z `~`). Na chwile obecną należy pamiętać, że formuły wyglądają mniej więcej tak: `Y~X`. 

```{r, fig.height=5, fig.width=15}
# Ustawiamy parametry graficzne R tak, aby narysować trójpanelowy wykres
par(mfrow = c(1,3))
# Tworzymy pierwszy wykres za pomocą składni formuły
plot(data$FEELBAD ~ data$ART)
# Przypisujemy model do zmiennej `model`
model <- lm(data$FEELBAD ~ data$ART)
# Model można przekazać jako argument dla funkcji `abline` i ona będzie wiedziała co i gdzie narysować
abline(model)

# Pozostałę dwa przypadki analogicznie
plot(data$FEELBAD ~ data$IDEAL)
abline(lm(data$FEELBAD ~ data$IDEAL))

plot(data$FEELBAD ~ data$SOCIAL)
abline(lm(data$FEELBAD ~ data$SOCIAL))
```

## Linia regresji

Równanie prostej dopasowywanej do danych w regresji liniowej ma postać:

$$\hat{Y} = bX + a$$

gdzie:

$\hat{Y}$ - przewidywana wartość Y

$b$ - współczynnik regresji (*slope* - współczynnik kierunkowy) 

$a$ - wyraz wolny (*intercept*)

$X$ - wartość zmiennej predyktora

```{r}
plot(FEELBAD ~ PEOPLE, data = data)
abline(lm(FEELBAD ~ PEOPLE, data = data))
summary(lm(FEELBAD ~ PEOPLE, data = data))
```

### Wyraz wolny (*intercept*)

* **definicja**: wartość $\hat{Y}$ kiedy $X$ przyjmuje wartość 0
* jego interpretacja zależy od tego, czy $X = 0$ ma jakąkolwiek sensowną interpretacje
* zwykle nie ma sensownej interpretacji i ma tylko tę matematyczną (ogromna i niepraktyczna ekstrapolacja z naszych danych - pomyśl o wadze 0kg!)
* zawsze można **wycentrować** nasz predyktor wokół średniej - wtedy uzyskujemy sensowną interpretację - $\hat{Y}$ dla wartości oczekiwanej (czyli średniej!) $X$
* **wycentrowanie** nie ma żadnych skutków dla współczynnika kierunkowego oraz współczynnika korelacji

### Współczynnik kierunkowy (*slope*)

* zmiana $\hat{Y}$ związana ze zmianą $X$ o jedną jednostkę
* definicja ta mówi nam że ma on sensowną interpretację (np. jeżeli mówimy o regresji dochodu z liczby lat edukacji, to współczynnik kierunkowy powie nam jaka różnica w dochodzie jest związana z każdym dodatkowym rokiem edukacji)

### Standaryzowany współczynnik regresji

* to co współczynnik kierunkowy, tylko że obie zmienne są wystandaryzowane
* możemy policzyć za pomocą `lm.beta` z pakietu `QuantPsyc`, ale równie dobrze możemy zrobić to sami mnożąc przez iloraz wariancji

### Korelacja a standaryzowany współczynnik regresji

* jeżeli mamy jeden predyktor (tak jak w tych przykładach, których się zajmujemy) to jest to ta sama wartość

**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ółczynnikiem 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 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

Kontunuując prace z danymi Covy i Giffina możemy spytać, czy istnieje zależność między odpowiedziami na pytania o normatywną wymiar doświadczenia *guilty pleasure* ("I feel like i should not enjoy this work" oraz "I feel like my ojoyment of this work is wrong and misguided") a tym, jak źle czują się z tym badani.

```{r}
data$NORM <- apply(data[,c("SHOULDNOT", "MISGUIDED")], 1, mean)
```

Nastepnie dopasujmy model liniowy do naszych danych.

```{r}
model <- lm(FEELBAD ~ NORM, data= data)
summary(model)
```

Zobaczmy jak wygląda wykres punktowy (wykres rozrzutu) dla zmiennej objaśniającej i objaśnianej.

```{r}
plot(FEELBAD ~ NORM, data = data)
abline(model)
```

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(model)
```