--- title: Korelacja. Regresja format: html: embed-resources: true code-copy: true lang: pl theme: cosmo format-links: false pdf: pdf-engine: xelatex include-in-header: text: | \usepackage{fvextra} \DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\},breaklines,breakanywhere} \RecustomVerbatimEnvironment{verbatim}{Verbatim}{breaklines,breakanywhere} lang: pl engine: knitr editor: source execute: error: false --- W tej części kursu będziemy zajmować się korelacją oraz regresją. Na początek kilka uwag wstępnych. W statystyce zasadniczo odróżnia się te dwa pojęcia na poziomie konceptualnym. Rozróżnienie to nie jest konsekwentnie przez wszystkich stosowane, warto jednak je znać (dla spokoju ducha!). **Regresja** - metoda pozwalająca na zbadanie związku pomiędzy zmiennymi i wykorzystanie tej wiedzy do przewidywania nieznanych wartości jednych wielkości na podstawie znajomości wartości innych. Poszukuje się zwązku między jedną (lub więcej) zmienną objaśniającą lub niezależną $X$ a zmienną objaśnianą lub zależną $Y$. W **regresji** zmienna $X$ (objaśniająca) jest w pełni kontrolowana przez eksperymentatora i pozbawiona elementu losowości. W **korelacji** obie zmienne sa zmiennymi losowymi. W praktyce, tak jak wspomnieliśmy, różnica ta jest trochę zatarta, ale warto o niej pamiętać. # Wykres punktowy (*scatterplot*) Na początek przyjrzyjmy się sposobowi wizualizacji relacji między dwiema zmiennymi. Nie zdziwi Państwa informacja, że znany i (przynajmniej przez niektórych) lubiany wykres punktowy (zwany takze wykresem rozrzutu) świetnie nadaje się do ilustrowania tego typu danych. Przypomnijmy więc sobie, w jaki sposób w R możemy narysować wykres punktowy. ## Podstawy Domyślnie, jeżeli wywołamy funkcje `plot` i przekażemy jej jako argumenty dwa wektory typu `numeric`, to R narysuje wykres punktowy. W tym przypadku reprodukujemy wykresy znajdujące się w rozdziale 9 podręcznika Howella. Każdy z punktów na wykresie reprezentuje jeden kraj. Pierwszy wykres przedstawia relację między liczbą lekarzy (zmienna objaśniająca) a śmiertelnością noworodków (zmienna objaśniana), drugi między wydatkami na służbę zdrowia (zmienna objaśniająca) a oczekiwaną długością życia (zmienna objaśniana), trzeci zaś między promieniowaniem słonecznym (zmienna objaśniająca) a zachorowaniem na nowotwory (zmienna objaśniana). ```{r fig.width=7.5, fig.height=4} par(mfrow = c(1,3)) d1 <- read.table('Fig9-1a.dat', header = TRUE) plot(d1$Physicians, d1$InfMort) d2 <- read.table('Fig9-1b.dat', header = TRUE) plot(d2$Expend, d2$LifeExp) d3 <- read.table('Fig9-1c.dat', header = TRUE) plot(d3$Radiation, d3$Cancer) ``` ## Dlaczego zawsze należy wizualizować dane? Datasaurus Datasaurus (czyli dinozaur ułożony z punków stworzony przez badaczy z firmy Autocad; więcej informacji na: https://www.autodeskresearch.com/publications/samestats) jest dowodem na to, że zawsze powinniśmy wizualizować nasze dane przed przeprowadzeniem analiz statystycznych. Poniżej znajduje się ilustracja 13 rozkładów, które mają takie same: - średnie X - średnie Y - odchylenie standardowe X - odchylenie standardowe Y - współczynnik korelacji między X a Y - współczynnik kowariancji między X a Y A mimo to są diametralnie różnymi rozkładami! Gdybyśmy patrzyli wyłącznie na gołe statystyki liczbowe, to nigdy byśmy nie zobaczyli, że te zbiory danych są diametralnie różne. Oto wykresy: ```{r fig.width=15, fig.height=15} options(repr.plot.width=15, repr.plot.height=15) library(ggplot2) ggplot(data = datasauRus::datasaurus_dozen) + geom_point(aes(x = x, y = y)) + facet_wrap(~dataset, ncol=4) ``` A tutaj znajdą Państwo statystyki deskryptywne! ```{r message=F, warning=F} library(dplyr) group_by(datasauRus::datasaurus_dozen, dataset) %>% summarise( "Korelacja" = cor(x, y), "Kowariancja" = cov(x, y), "Średnia X" = mean(x), "Średnia Y" = mean(y), "Odchylenie standardowe X" = sd(x), "Odchylenie standardowe Y" = sd(y), ) ``` # 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 informacje o wywołaniu funkcji. Jeżeli chcemy się dowiedzieć czegoś więcej, musimy obiekt ten przekazać funkcji `summary`. ## Wywołanie funkcji `lm` Wywołanie funkcji `lm` na danych zwraca nam obiekt klasy `lm`. Możemy fo przypisać do zmiennej (często spotykaną w R konwencją jest nazwa `fit` dla modelu) albo po prostu wyświetlić. ```{r} # Wywołanie funkcji `lm` zwraca nam odpowiednio dopasowany model. lm(d1$Physicians ~ d1$InfMort) ``` ## Wywołanie funkcji `summary` na obiekcie zwracanym przez `lm` Tak, jak powiedzieliśmy sobie wczesniej, funkcja `summary` wywołana na obiekcie zwróconym przez funkcję `lm` pozwala nam się dowiedzieć więcej o stworzonym przez nas modelu. W szczególności zaś pozwala nam poznać różne jego parametry oraz przeprowadza szereg testów statystycznych (na istotność poszczególnych współczynników w regresji oraz na dopasowanie całego modelu). ```{r} summary(lm(d1$Physicians ~ d1$InfMort)) ``` ### Dodawanie linii regresji Moglibyśmy ręcznie ,,wydobyć'' wartości niezbędne do dodania linii regresji (są w końcu drukowane przez R!), ale możemy skorzystać z faktu, że nasz model stworzony za pomocą funkcji `lm` możemy bezpośrednio przekazać jako argument do funkcji `abline`. Funkcja ta zaś automatycznie dorysuje do istniejącego wykresu linię regresji. 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} # Ustawiamy parametry graficzne R tak, aby narysować trójpanelowy wykres par(mfrow = c(1,3)) # Tworzymy pierwszy wykres za pomocą składni formuły (tzn. składni z tyldą: Y ~ X) plot(d1$Physicians ~ d1$InfMort) # Przypisujemy model do zmiennej `fit` fit <- lm(d1$Physicians ~ d1$InfMort) # Model można przekazać jako argument dla funkcji `abline` i ona będzie wiedziała, co i gdzie narysować abline(fit) # Pozostałę dwa przypadki analogicznie plot(d2$LifeExp ~ d2$Expend) abline(lm(d2$LifeExp ~ d2$Expend)) plot(d3$Radiation ~ d3$Cancer) abline(lm(d3$Radiation ~ d3$Cancer)) ``` # Kowariancja i współczynnik korelacji *r* Pearsona ## Kowariancja Wzór na współczynnik kowariancji między dwiema zmiennymi wygląda następująco: $$ Cov(X,Y) = E[(X-E(X)) (Y-E(Y))] $$ * bardzo podobna do wariancji (gdyby za $Y$ podstawić $X$ byłby to wzór na wariancje) * miara współzmienności dwóch zmiennych * jeżeli dużym $X$ (w sensie odległości od średniej) towarzyszą duże $Y$, to kowariancja będzie wysoka i dodatnia, jeśli małe $Y$ to będzie ujemna, jeśli raz takie, a raz takie (brak korelacji) to będą się znosić a kowariancja będzie wynosić około 0 ## Korelacja Problem z kowariancją jest taki, że jej wielkość zależy od tego, jakie jednostki mają nasze zmienne oraz w szczególności od tego, jakie mają odchylenie standardowe. Może to utrudnić porównania i interpretacje takiej wartości. Obliczenie współczynnika korelacji jest sposobem na poradzenie sobie z tym problem. Można o nim myśleć jako o "wystandaryzowanym" współczynniku kowiariancji. Wzór na korelację w populacji ($\rho$) wygląda tak: $$ \rho = \frac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}} $$ * przyjmuje wartości z przedziału $[-1, 1]$ * jest miarą liniowej zależności między zmiennymi losowymi $X$ i $Y$ * jeżeli dysponujemy próbą, to możemy wyliczyć współczynnik korelacji Pearsona ($r$): $$ r = \frac{\sum_{i=1}^{n} (X_i - \bar{X})(Y_i - \bar{Y})}{\sqrt{\sum_{i=1}^{n}(X_i - \bar{X})^2 (Y_i - \bar{Y})^2 }}$$ albo prościej: $$r = \frac{cov_{XY}}{s_x s_y}$$ * na podstawie jego wartości możemy ocenić siłę związku prostoliniowego między cechami $X$ a $Y$: Spróbujmy sprawdzić teraz, czy R oblicza oba współczynniki zgodnie ze wzorami! Zaczniemy od kowariancji: ```{r} print('Kowariancja') sum((d1$Physicians - mean(d1$Physicians)) * (d1$InfMort - mean(d1$InfMort)))/ (nrow(d1)-1) # wzór kowariancji wpisany ręcznie cov(d1$Physicians, d1$InfMort) # funkcja wbudowana w R ``` A teraz sprawdzimy to samo dla korelacji: ```{r} print('Korelacja') cor(d1$Physicians, d1$InfMort) # funkcja wbudowana w R (sum((d1$Physicians - mean(d1$Physicians))*(d1$InfMort - mean(d1$InfMort)))/ (nrow(d1)-1))/(sd(d1$Physicians)*sd(d1$InfMort)) ``` ## Interpretacja wartości *r* Interpretując współczynnik korelacji Pearsona musimy pamiętać, że właściwa interpretacja będzie zależeć od tego, w jakiej dziedzinie się aktualnie znajdujemy. Dla fizyka inny współczynnik korelacji będzie uznany za "duży" niż dla psychologa społecznego. Poniżej znajdują się jednak pewne "zasady kciuka" dotyczące interpretacji współczynnika $r$ w naukach psychologicznych: - $r = 0$ - współzależność nie występuje, brak korelacji, zmienne są nieskorelowane - $0 < \mid r \mid < 0.3$ - słaby stopień współzależności - $0.3 \leq \mid r \mid < 0.5$ - średni stopień współzależności - $0.5 \leq \mid r \mid < 0.7$ - znaczny stopień współzależności - $0.7 \leq \mid{r}\mid < 0.9$ - wysoki stopień współzależności - $\mid r \mid \geq 0.9$ - bardzo wysoki stopień współzależności - $\mid r \mid = 1$ - współzależność całkowita (ścisłość) ## Macierze korelacyjne - wizualizacje Czasami mamy do czynienia z więcej niż dwiema zmiennymi. W takim przypadkach zdarza się, że chcielibyśmy zbadać korelacje między wszystkimi kombinacjami dwóch zmiennych jednoczesnie. W tym celu możemy stworzyć macierz korelacyjną (tak samo tworzy się macierz kowariancji - przyda się Państwu na Statystyce II przy PCA). Tutaj zrobimy to (może nieco nudno) na losowych danych pochodzacych z rozkładu normalnego. ```{r} # Wylosujemy 350 obserwacji z rozkładu normalnego a następnie podzielimy je na 10 kolumn # W ten sposób zasymulujemy sytuacje zbioru danych z 10 zmiennymi m <- matrix(rnorm(350), 35, 10) # Przekształcamy macierz w ramkę danych m <- as.data.frame(m) # Tworzymy macierz korelacyjną, dla kowariancji będzie to `cov` cor(m) ``` Jednym ze sposobów wizualizacji takich danych jest "mapa cieplna", w której intensywność koloru (albo jakiś jego inny parametr) oznacza siłę korelacji (tutaj - im jaśniejszy tym wyższy współczynnik korelacji, im ciemniejszy tym niższy). ```{r} library('viridis') heatmap(cor(m), symm = TRUE, Rowv = NA, col=viridis(256)) ``` Wersja podstawowa produkowana przez funkcję `heatmap` wygląda w tym przypadku tak. ```{r} heatmap(cor(m)) # wersja z klastrami, ale one nam są niepotrzebne ``` Pewną ciekawą odmianą mapy cieplnej jest wykres z elipsami, na którym widać siłę (jak bardzo spłaszczona jest elipsa) oraz kierunek (w którą stronę ,,kopnięta'' jest elipsa) korelacji między zmiennymi. ```{r} library(ellipse) # `type` mówi o tym, czy chcemy umieszczać na wykresie całą macierz czy tylko jej część # `diag` mówi, czy rysować przekątną (na przekątnej r = 1, więc nie jest za bardzo ciekawa...) plotcorr(cor(m), type = 'lower', diag = TRUE) ``` W przypadku wylosowanych obserwacji może być trudno zrozumieć, jak działać mają rysowane prez `plotcorr` elipsy. Być może lepiej widać to na znanym Państwu już zbiorze dotyczącym irysów (`iris`). ```{r} plotcorr(cor(iris[,1:4]), type = 'lower', diag = TRUE) ``` Ostatnim sposobem generowania ładnie wyglądających wizualizacji korelacji między wieloma zmiennymi, który omówimy, są funkcje `corrplot` i `corrplot.mixed` z pakietu (!) `corrplot`. Funkcje ta ma bardzo wiele możliwości dostosowania wyglądu wykresu - zachęcam do samodzielnej eksploracji! ```{r} library(corrplot) corrplot.mixed(cor(iris[,1:4]), lower = "ellipse", upper = "number") ``` # Gra - zgadnij korelacje Jeśli będziemy w przyszłości pracować z danymi, to warto nabyć pewną intuicję dotyczącą tego, jaki wizualny "rozrzut" punktów odpowiada jakiemu współczynnikowi korelacji. Dobrym ćwiczeniem jest dostępna tutaj gra, w której możemy przetestować i "poprawić" swoją intuicję. http://guessthecorrelation.com/ Możemy podobną grę przeprowadzić za pomocą R. Dobrym punktem wyjścia jest funkcja `rmvnorm` z pakietu `mvtnorm`, która pozwala nam losować skorelowane ze sobą zmienne. ```{r} library(mvtnorm) sym <- rmvnorm(100, mean = c(1,1), sigma = matrix(c(1,1,1,1),2,2)) plot(sym) cov(sym) cor(sym) ``` ```{r} sym <- rmvnorm(100, mean = c(1,1), sigma = matrix(c(3,2,2,3),2,2)) plot(sym) cov(sym) cor(sym) ``` # Linia regresji Na koniec powróćmy do regresji liniowej. Jak powszechnie wiadomo, 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 ## Zadanie Wagner, Compas i Howell (1988) badali związek między stresem a zdrowiem psychicznym wśród studentów pierwszego roku koledżu. Używając opracowanego przez siebie narzędzia do mierzenia częstotliwości, odczuwanej wagi i pożądaności niedawnych wydarzeń życiowych, stworzyli miarę negatywnych zdarzeń w życiu. Miara ta służyła do pomiaru środowiskowego i społecznego stresu odczuwanego przez badanych. Oprócz tego poprosili swoich badanych (studentów), aby wypełnili *Hopkins Symptom Checklist*, która to służy do pomiaru występowania lub brak występowania 57 psychologicznych symptomów zaburzeń zdrowia psychicznego. Rozpoczniemy od wczytania oraz wizualizacji danych. Za pomocą funkcji `abline` dodamy do wykresu punktowego linię regresji. ```{r} df <- read.table('Tab9-2.dat', header = TRUE) head(df) plot(lnSymptoms ~ Stress, data = df) abline(lm(lnSymptoms ~ Stress, data = df)) ``` Następnie przyjrzymy się bliżej stworzonemu przez nas modelowi. ```{r} summary(lm(lnSymptoms ~ Stress, data = df)) ``` (Dygresja: w pakiecie `rms` znajduje się funkcja do przeprowadzania regresji liniowej zwracająca nieco inne informacje niz ta domyślnie wbudowana w R. Można ją przetestować po wczytaniu pakietu `rms` (`library(rms)`) wywołując polecenie `ols(lnSymptoms ~ Stress, data = df)`) ## 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 $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ść # Testowanie hipotez dla współczynnika korelacji $r$ Pearsona Najczęściej będziemy sprawdzać czy zmienne $X$ i $Y$ są skorelowane. W takim wypadku nasze hipotezy będą przestawiać się w następujący sposób. Zaczniemy od hipotezy alternatywnej. $$ H_A: \rho \neq 0 $$ To znaczy, że hipoteza alternatywna głosi, że korelacja w populacji ($\rho$) jest różna od zera. Hipoteza zerowa (tę będziemy obalać!) głosi więc, że korelacja w populacji wynosi 0: $$ H_0: \rho = 0 $$ Statystyka testowa wygląda w takim przypadku tak: $$T = \frac{r}{\sqrt{1-r^2}}\sqrt{n-2}$$ która ma przy prawdziwości $H_0$ rozkład $t$ o n-2 stopniach swobody. Możemy to z łatwością sprawdzić za pomocą prostego eksperymentu symulacyjnego. ```{r} n <- 35 t_stat = replicate(10000, { r = cor(rnorm(n), rnorm(n)) (r/sqrt(1-r^2))*sqrt(n-2)}) hist(t_stat, freq = FALSE) curve(dt(x,n-2), add = TRUE) ``` Zróbmy proste ćwiczenie. Na początek stwórzmy macierz z losowymi liczbami pochodzącymi z rozkładu normalnego (10 kolumn po 35 liczb) a następnie obliczmy statystykę testową $t$ dla korelacji między nimi. ```{r} # Generujemy zbiór danych m <- matrix(rnorm(350), 35, 10) # Tworzymy macierz korealcji c_m <- cor(m) # Tworzymy macierz statystyk t t_m <- (c_m/sqrt(1-c_m^2)) * (sqrt(35-2)) t_m ``` Mając taką macierz jesteśmy w stanie dla każdej komórki (= dla każdej kombinacji dwóch zmiennych) ocenić, czy odrzucamy hipotezę zerową czy nie. Wystarczy posłużyć się odpowiednim rozkładem $t$ o $n-2$ stopniach swobody. Za pomocą funkcji poniżej możemy łatwo stwierdzic, które hipotezy zerowe odrzucamy. ```{r} t_m <= qt(0.025, 33) | t_m >= qt(0.975, 33) ``` A nawet policzyć, ile razy odrzuciliśmy hipotezę zerową! ```{r} print('Liczba statystycznie istotnych korelacji między wygenerowanymi próbami') sum(t_m <= qt(0.025, 33) | t_m >= qt(0.975, 33)) - ncol(t_m) ``` (Uwaga! Test statystyczny dla współczynnika korelacji możemy również wykonać za pomocą wbudowanej w R funkcji `cor.test`. Funkcja ta jednak przyjmuje tylko dwa wektory! ```{r} cor.test(m[, 1], m[, 2]) ``` ) Poniżej widzą Państwo wykres ilustrujący nasz eksperyment z losowaniem. Pod wykresem wypisujemy liczbę istotnych statystycznie współczynników korelacji poza przekątną (czyli pomijamy korelacje zmiennych z samymi sobą). Wokół jakiej liczby oscylowałaby ta wartość, gdybyśmy powtarzali nasz eksperyment? Dlaczego? ```{r} m <- matrix(rnorm(350), 35, 10) c_m <- cor(m) t_m <- (c_m/sqrt(1-c_m^2)) * (sqrt(35-2)) reject <- (t_m <= qt(0.025, 33) | t_m >= qt(0.975, 33)) heatmap(matrix(as.numeric(reject), 10,10), Rowv = NA, Colv = NA, col = viridis(2)) print(sum(reject) - 10) # 10 leży na przekątnej ```