library(ggplot2)
library(GGally)
library(Hmisc)
library(psych)
library(PerformanceAnalytics)
library(car)
library(corrplot)

1 Podstawowe informacje

  • Korelacja oznacza związek pomiędzy zmiennymi.

  • Analiza korelacji służy do “wychwycenia” czy zachodzi związek pomiędzy dwiema (lub więcej) zmiennymi.

  • Miarą korelacji jest współczynnik korelacji

  • Współczynnik korelacji dostarcza informacji o tym jaka jest siła związku (wartość współczynnika) oraz jaki jest kierunek związku (znak).

  • Dla każdego współczynnika korelacji należy także obliczyć jego istotność statystyczną, stosujac jeden z testów istotności przeznaczonych dla współczynników korelacji.

    • Hipoteza zerowa: ρ x,y=0
    • Hipoteza alternatywna: ρ x,y≠0 lub ρ x,y<0 lub ρ x,y>0

  • Zależność między zmiennymi może mieć charakter liniowy lub krzywoliniowy.

2 Korelacja a przyczynowość

  • Korelacja nie wskazuje na istnienie związku przyczynowo-skutkowego pomiędzy zmiennymi.

  • Innymi słowy: Istnienie korelacji liczbowej nie potwierdza, że jedno zjawisko powoduje drugie.

    • A może powodować B
    • B może powodować A
    • A lub B może być wywołane przez C
    • Zależność między A i B może być przypadkowa.

  • http://www.tylervigen.com/spurious-correlations

3 Kilka ważnych informacji

  • Najważniejsza jest isotność korelacji. Niepotrzebna nam korelacja nawet bardzo wysoka, jeśli nie jest istotna statystycznie.
  • Wartość współczynnika nawet bliska 0 nie zawsze oznacza brak zależności. Może oznaczać jedynie brak zależności liniowej.
  • Wielkość współczynnika podlega wpływom wartości skrajnych i odstających.

4 Korelacja liniowa

  • Miarą korelacji liniowej jest współczynnik korelacji Pearsona.

  • Współczynniki korelacji przyjmują wartości z przedziału od -1,00 do +1,00.

    • Wartość -1,00 - reprezentuje doskonałą korelację ujemną (współzależność pomiędzy zmiennymi kształtująca się w taki sposób, że gdy wartości jednej zmiennej wykazują tendencję rosnącą, wówczas wartości drugiej zmiennej wykazują tendencję malejącą)

    • wartość +1,00 - reprezentuje doskonałą korelacją dodatnią (współzależność pomiędzy zmiennymi przedstawia się w taki sposób, że gdy wartości jednej zmiennej wykazują tendencję wzrastającą, wówczas wartości drugiej zmiennej także wykazują tendencję wzrastającą).

    • Wartość 0.00 wyraża brak korelacji.

4.1 Jak silna jest korelacja?

Do opisu i interpretacji istotnej korelacji pomocne może być przyjęcie pewnej skali określającą siłę związku. Nie ma jednej przyjętej skali. Poniżej przedstawiam jedną z nich:

  • 0 - brak korelacji
  • 0,1 do 0,3 - słaba korelacja
  • 0,3 do 0,5 - przeciętna korelacja
  • 0,5 do 0,7 - korelacja wysoka
  • 0,7 do 0,9 - korelacja bardzo wysoka
  • 0,9 do 1 - korelacja prawie pełna

4.2 Demonstracja dla współczynnika korelacji liniowej

install.packages("TeachingDemos")
library("TeachingDemos")
put.points.demo(x = NULL, y = NULL, lsline = TRUE)

#Używając opcji Add Point dodaj punkty w oknie wykresu
#Zwróć uwagę jak zmienia się wartość współczynnika korelacji (r).

Zgadnij wartość współczynnika korelacji - https://gallery.shinyapps.io/correlation_game/

5 Testy korelacji

  • test korelacji liniowej Pearsona

    • stosowany gdy zmienne mają zależnośc liniową
    • zmienne mają rozkład normalny

  • test korelacji rang Spearman

    • stosowany gdy naruszone jest założenie o normalności rozkładu
      (np. gdy istnieją wartości odstające)

6 PRZYKŁAD: Czy istnieją zależności liniowe pomiędzy zmiennymi w zbiorze cardata?

cardata <- read.csv("dane/cardata2.csv")

6.1 Określenie zależności pomiędzy dwiema zmiennymi

6.1.1 Wykres rozrzutu

ggplot(cardata, aes(x = mpg, y = horsepower)) + geom_point() + theme_bw()
## Warning: Removed 5 rows containing missing values (geom_point).

#Z linią najlepszego dopasowania
ggplot(cardata, aes(x = mpg, y = horsepower)) + geom_point() +geom_smooth()

#Z linią najlepszego dopasowania dla zależności liniowej
ggplot(cardata, aes(x = mpg, y = horsepower)) + geom_point() +geom_smooth(method = 'lm')

6.1.2 Sprawdzenie zgodności z rozkładem normalnym

Czy zmienna horsepower oraz mpg mają rozkład normalny lub zbliżony do normalnego?

p1 <- ggplot(cardata, aes(x = mpg)) + geom_histogram() + theme_bw()
p2 <- ggplot(cardata, aes(x = horsepower)) + geom_histogram() + theme_bw()
gridExtra::grid.arrange(p1, p2, ncol=2)

shapiro.test(cardata$mpg)
## 
##  Shapiro-Wilk normality test
## 
## data:  cardata$mpg
## W = 0.95549, p-value = 0.001034
shapiro.test(cardata$horsepower)
## 
##  Shapiro-Wilk normality test
## 
## data:  cardata$horsepower
## W = 0.94986, p-value = 0.0005026

6.1.3 Współczynnik korelacji liniowej Pearsona

cor(cardata$mpg, y = cardata$horsepower, use = "complete.obs", method = c("pearson"))
## [1] -0.7972331

6.1.4 Test istotności dla współczynnika korelacji liniowej Pearsona

cor.test(cardata$mpg, y = cardata$horsepower, method = c("pearson"))
## 
##  Pearson's product-moment correlation
## 
## data:  cardata$mpg and cardata$horsepower
## t = -13.468, df = 104, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8575719 -0.7152504
## sample estimates:
##        cor 
## -0.7972331

6.1.5 Współczynnik korelacji rang Spearmana

cor(cardata$mpg, y = cardata$horsepower, use = "complete.obs", method = c("spearman"))
## [1] -0.8342363

6.1.6 Test istotności dla współczynnika korelacji rang Spearmana

cor.test(cardata$mpg, y = cardata$horsepower, method = c("spearman"))
## Warning in cor.test.default(cardata$mpg, y = cardata$horsepower, method =
## c("spearman")): Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  cardata$mpg and cardata$horsepower
## S = 364068, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.8342363

Czy między zmiennymi horsepower oraz mpg istnieje zależność liniowa? Która z metod (współczynnik korelacji Pearsona czy współczynnik korelacji rang Spearmana) powinna zostać użyta do określenia zależności między zmiennymi horsepower oraz mpg? Dlaczego?

Wykonaj wykres rozrzutu, oblicz współczynnik korelacji oraz sprawdź czy wynik jest istotny statystycznie dla zmiennych mpg oraz accell. Która miara korelacji (współczynnik korelacji Pearsona czy współczynnik korelacji rang Spearmana) powinna zostać użyta w tym przykładzie?

#wykres rozrzutu
ggplot(cardata, aes(x = mpg, y = accel)) + geom_point()

#histogramy zmiennych
p1 <- ggplot(cardata, aes(x = mpg)) + geom_histogram() + theme_bw()
p2 <- ggplot(cardata, aes(x = accel)) + geom_histogram() + theme_bw()
gridExtra::grid.arrange(p1, p2, ncol=2)

#testowanie zgodnosci z rozkladem normalnym 
shapiro.test(cardata$mpg)
shapiro.test(cardata$accel)

#korelacja
cor(cardata$mpg, y = cardata$accel, use = "complete.obs", method = c("spearman"))
cor.test(cardata$mpg, y = cardata$accel, method = c("spearman"))

6.2 Określanie korelacji dla wielu zmiennych

6.2.1 Macierz wykresów korelacji

pairs(~mpg+displace+accel+horsepower+weight, data=cardata, main="Simple Scatterplot Matrix")

Interpretacja macierzy wykresów rozrzutu

Wykres w lewym dolnym rogu pokazuje zależność między mpg (oś x) a displace (oś y).

6.2.2 Macierz współczynników korelacji - Współczynnik korelacji liniowej Pearsona

cardata2 <- cardata[,c(-6,-8,-9,-10)]

Funkcja cor() pozwala na obliczenie macierzy współczynników korelacji. Podaje ona wartość współczynnika, ale nie wskazuje, czy wynik jest istotny statystycznie.

Funkcja rcorr() z pakietu Hmisc wyświetla zarówno macierz współczynników korelacji, jak i wartość p wskazującą czy wynik jest istotny statystycznie.

korelacja <- cor(cardata2, y = NULL, use = "complete.obs", method = c("pearson"))
round(korelacja,2)
##              mpg cylinders displace horsepower accel weight price
## mpg         1.00     -0.69    -0.74      -0.80  0.25  -0.82  0.08
## cylinders  -0.69      1.00     0.94       0.78 -0.23   0.83 -0.02
## displace   -0.74      0.94     1.00       0.82 -0.23   0.90  0.03
## horsepower -0.80      0.78     0.82       1.00 -0.48   0.79 -0.06
## accel       0.25     -0.23    -0.23      -0.48  1.00  -0.01  0.26
## weight     -0.82      0.83     0.90       0.79 -0.01   1.00  0.17
## price       0.08     -0.02     0.03      -0.06  0.26   0.17  1.00

W poniższym przykładzie pierwsza macierz zawiera współczynnik korelacji, druga liczbę obiektów a trzecia wartość poziomu istotności p. Wartość jest istotna statystycznie jeśli p jest mniejsze od założonego poziomu isotntości (np. 0,05)

library(Hmisc)
kor <- rcorr(as.matrix(cardata2), type = "pearson")
kor
##              mpg cylinders displace horsepower accel weight price
## mpg         1.00     -0.68    -0.74      -0.80  0.24  -0.83  0.06
## cylinders  -0.68      1.00     0.93       0.78 -0.23   0.82 -0.03
## displace   -0.74      0.93     1.00       0.82 -0.22   0.90  0.03
## horsepower -0.80      0.78     0.82       1.00 -0.48   0.78 -0.04
## accel       0.24     -0.23    -0.22      -0.48  1.00  -0.01  0.25
## weight     -0.83      0.82     0.90       0.78 -0.01   1.00  0.18
## price       0.06     -0.03     0.03      -0.04  0.25   0.18  1.00
## 
## n
##            mpg cylinders displace horsepower accel weight price
## mpg        110       110      110        106   110    110   110
## cylinders  110       111      111        107   111    111   111
## displace   110       111      111        107   111    111   111
## horsepower 106       107      107        107   107    107   107
## accel      110       111      111        107   111    111   111
## weight     110       111      111        107   111    111   111
## price      110       111      111        107   111    111   111
## 
## P
##            mpg    cylinders displace horsepower accel  weight price 
## mpg               0.0000    0.0000   0.0000     0.0122 0.0000 0.5114
## cylinders  0.0000           0.0000   0.0000     0.0162 0.0000 0.7864
## displace   0.0000 0.0000             0.0000     0.0202 0.0000 0.7867
## horsepower 0.0000 0.0000    0.0000              0.0000 0.0000 0.6463
## accel      0.0122 0.0162    0.0202   0.0000            0.9488 0.0091
## weight     0.0000 0.0000    0.0000   0.0000     0.9488        0.0596
## price      0.5114 0.7864    0.7867   0.6463     0.0091 0.0596

W celu łatwiejszego przeglądu wyników, powyższe macierze można przekształcić do ramki danych z czterema kolumnami (col1 - nazwa zmiennej 1, col2 - nazwa zmiennej 2, col3 - współczynnik korelacji, col4 - wartość p).

# ++++++++++++++++++++++++++++
# flattenCorrMatrix
# ++++++++++++++++++++++++++++
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
}
kor_df <- flattenCorrMatrix(kor$r, kor$P)
head(kor_df)
##         row     column        cor            p
## 1       mpg  cylinders -0.6813499 2.220446e-16
## 2       mpg   displace -0.7426211 0.000000e+00
## 3 cylinders   displace  0.9346452 0.000000e+00
## 4       mpg horsepower -0.7972331 0.000000e+00
## 5 cylinders horsepower  0.7754365 0.000000e+00
## 6  displace horsepower  0.8169579 0.000000e+00

Z powyższej tabeli można wyselekcjonować, tylko te przypadki, dla których wartóść korelacji jest istotna

kor_df_sel <- kor_df[kor_df$p<0.05,]

6.2.3 Macierz współczynników korelacji - Współczynnik korelacji rang Spearmana

korelacja_rang <- cor(cardata2, y = NULL, use = "complete.obs", method = c("spearman"))
round(korelacja_rang,2)
##              mpg cylinders displace horsepower accel weight price
## mpg         1.00     -0.76    -0.84      -0.83  0.13  -0.85  0.09
## cylinders  -0.76      1.00     0.87       0.74 -0.17   0.83 -0.02
## displace   -0.84      0.87     1.00       0.85 -0.15   0.92  0.03
## horsepower -0.83      0.74     0.85       1.00 -0.35   0.77  0.00
## accel       0.13     -0.17    -0.15      -0.35  1.00   0.07  0.16
## weight     -0.85      0.83     0.92       0.77  0.07   1.00  0.10
## price       0.09     -0.02     0.03       0.00  0.16   0.10  1.00
kor <- rcorr(as.matrix(cardata2), type = "spearman")

6.2.4 Zależność pomiedzy zmiennymi w zbiorze cardata:

Interpretacja wykresów rozrzutu – dane cardata

  • horsepower-weight (moc silnika a waga samochodu) – zależność wprost proporcjonalna (im większa moc silnika, tym cięższy silnik a zatem i większa waga samochodu)

  • horsepower-mpg (moc silnika a ilość mil przejechanych na jednym galonie)– zależność odwrotnie proporcjonalna – wraz ze wzrostem mocy silnika maleje ilość mil jaką możemy przejechać na jednym galonie tj. im większa moc silnika tym większe zużycie paliwa.

  • displace-mpg (pojemność silnika a ilość mil przejechanych na jednym galonie paliwa)- zależność krzywoliniowa

  • weight-mpg (waga a ilość mil przejechanych na jednym galonie) – obserwujemy wyraźne odstępstwa od zależności liniowej.

7 Porównanie wyników dla współczynnika korelacji liniowej oraz współczynnika rang Spearmana

  • horsepower-weight (moc silnika a waga samochodu) – Zależność liniowa wprost proporcjonalna
##  Pearsona  Spermana 
## 0.7837807 0.7686775
  • horsepower-mpg (moc silnika a ilość mil przejechanych na jednym galonie)– zależność odwrotnie proporcjonalna
##   Pearsona   Spermana 
## -0.7972331 -0.8342363
  • displace-mpg (pojemność silnika a ilość mil przejechanych na jednym galonie paliwa)- zależność krzywoliniowa
##   Pearsona   Spermana 
## -0.7426211 -0.8401613
  • weight-mpg (waga a ilość mil przejechanych na jednym galonie) – obserwujemy wyraźne odstępstwa od zależności liniowej.
##   Pearsona   Spermana 
## -0.8252383 -0.8562282

Współczynnik korelacji liniowej Pearsona jest miarą zależności liniowej. Jeśli zależność nie ma charakteru liniowego współczynnik ten będzie bardzo słabą miarą zależności. W takim wypadku stosuje się inną miarę – współczynnik korelacji rangowej Spearmana. Jeśli wartość współczynnika korelacji rangowej jest wyraźnie wyższa od wartości współczynnika korelacji liniowej tzn. że mamy do czynienia z zależnością krzywoliniową.

8 Wizualizacja korelacji

8.1 Funckja scatterplotMatrix() z pakietu car

library(car)
scatterplotMatrix(~mpg+displace+accel+horsepower+weight, data=cardata)

#Macierz korelacji w podziale na grupy
scatterplotMatrix(~mpg+displace+accel+horsepower+weight|origin, data=cardata)

8.2 Funkcja chart.Correlation() z pakietu PerformanceAnalytics

chart.Correlation(cardata2, histogram=TRUE, pch=19, method = "pearson")

chart.Correlation(cardata2, histogram=TRUE, pch=19, method = "spearman")

8.3 Funkcja pairs.panels() z pakietu psych

library(psych)
pairs.panels(cardata2, scale=TRUE)

8.4 Funkcja corrplot() z pakietu corrplot

library("corrplot")
kor <- rcorr(as.matrix(cardata2), type = "pearson")
corrplot(kor$r)

corrplot(kor$r, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

corrplot.mixed(kor$r, order="hclust", tl.col="black")

library(RColorBrewer)
corrplot(kor$r, type="upper", order="hclust", col=brewer.pal(n=8, name="RdYlBu"))

8.5 Funkcja ggcorr() i ggpairs() z pakietu GGally

library(GGally)
ggcorr(cardata2, nbreaks=8, palette='RdGy', label=TRUE, label_size=5, label_color='white')

library(GGally)
ggpairs(cardata2)

9 Zadanie do samodzielnego wykonania.

Używając danych z pakietu gapminder dla 2007 roku zwizualizuj korelację między zmiennymi pop, lifExp, gdpPercap. Oblicz korelację oraz wykonaj odpowiedni test. Którą miarę należy użyć w przypadku tych danych (współczynnik korelacji Pearsona czy współczynnik korelacji Spearmana)?