Analiza korelacji

Statystyczne i graficzne podstawy geowizualizacji

Author

Anna Dmowska, dmowska@amu.edu.pl 

library(tidyverse)
library(readxl)

1 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/

2 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)

3 Korelacja liniowa

set.seed(25)
x = rnorm(1000)
y = x + rnorm(1000)
df = data.frame(x, y)

Funkcja plot() pozwala na wykonanie prostego wykresu rozrzutu.

plot(y~x, df)

3.1 Zbadanie normalności rozkładu

Funkcja hist() pozwala na wykonie histogramów.

hist(df$x)

hist(df$y)

Funkcja shapiro.test() służy do wykonania testu Shapiro-Wilka sprawdzającego czy dane mają rozkład normalny. Wartość p-value > 0.05 wskazuje, że dane pochodzą z rozkładu normalnego.

shapiro.test(x)

    Shapiro-Wilk normality test

data:  x
W = 0.99811, p-value = 0.3298
shapiro.test(y)

    Shapiro-Wilk normality test

data:  y
W = 0.99854, p-value = 0.5832

Zmienna x oraz y mają rozkład normalny.

3.2 Współczynnik korelacji

Współczynnik korelacji obliczany jest z wykorzystaniem funkcji cor(). Argument method pozwala na określenie czy ma być obliczony współczynnik korelacji Pearsona (method = “pearson”), czy też współczynnik korelacji rang Spearmana (method = “Spearman”). Argument use = “complete.obs oznacza, że do obliczenia korelacji zostaną wybrane tylko te obserwacje, dla których jest podana wartość x oraz y.

cor(df$x, df$y,
    use = "complete.obs",
    method = "pearson")
[1] 0.679556
cor(df$x, df$y,
    use = "complete.obs",
    method = "spearman")
[1] 0.6636062

Która metoda korelacji powinna zostać użyta - korelacja liniowa Pearsona, czy korelacja rang Spearmana? Dlaczego?

3.3 Testy korelacji

Funkcja cor.test() wykonuje test korelacji. W wyniku testu:

  • p-value określa poziom istotności, jeśli p-value < 0.05 współczynnik korelacji jest istotny statystycznie
  • conf.int to 95% przedział ufności, w którym mieści się prawdziwa wartość współczynnika korelacji
  • sample estimates is to wartość współczynnika korelacji
cor.test(df$x, df$y,
    use = "complete.obs",
    method = "pearson")

    Pearson's product-moment correlation

data:  df$x and df$y
t = 29.263, df = 998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6447237 0.7115721
sample estimates:
     cor 
0.679556

Wynik testu korelacji wskazuje na istnieie istotnej korelacji między zmienną x oraz y.

4 Przykład. Analiza danych meteorologicznych.

Plik dane/meteo2020.xlsx zawiera miesięczne wartości dla 59 stacji synoptycznych dla 2020 roku:

  • id - id stacji
  • station - nazwa stacji
  • X, Y - współrzędne X, Y określające lokalizację stacji
  • wojewodztwo - nazwa województwa, w którym się stacja znajduje.
  • woj_id - TERYT województwa
  • mm - miesiąc
  • tmax_abs - miesięczna wartość temperatury maksymalnej [C]
  • tmax_mean - średnia miesięczna wartość temperatury maksymalnej [C]
  • tmin_abs - miesięczna wartość temperatury minimalnej [C]
  • tmin_mean - średnia miesięczna wartość temperatury minimalnej [C]
  • t2m_mean_mon - średnia temperatura miesieczna [C]
  • t5cm_min - minimalna temperatura przy gruncie [C]
  • rr_monthly - miesieczna suma opadow [mm]
  • rr_max_daily - maksymalna dobowa suma opadow [mm]
  • rain_days - liczba dni z opadem deszczu
  • snow_days - liczba dni z opadem śniegu
library(readxl)
meteo = read_excel("dane/meteo2020.xlsx")

4.1 Analiza korelacji między zmiennymi tmin_mean a tmax_mean

  • Wykres rozrzutu
plot(meteo$tmin_mean, meteo$tmax_mean)

  • Czy zmienne mają rozkład normalny?
shapiro.test(meteo$tmax_mean)

    Shapiro-Wilk normality test

data:  meteo$tmax_mean
W = 0.94879, p-value = 6.023e-15
shapiro.test(meteo$tmin_mean)

    Shapiro-Wilk normality test

data:  meteo$tmin_mean
W = 0.93303, p-value < 2.2e-16
  • wspolczynnik korelacji
cor(meteo$tmax_mean, meteo$tmin_mean,
    use = "complete.obs",
    method = "spearman")
[1] 0.8890966
  • test korelacji
cor.test(meteo$tmax_mean, meteo$tmin_mean,
    use = "complete.obs",
    method = "spearman")
Warning in cor.test.default(meteo$tmax_mean, meteo$tmin_mean, use =
"complete.obs", : Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  meteo$tmax_mean and meteo$tmin_mean
S = 6559830, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.8890966

5 Analiza korelacji dla wielu zmiennych

sel = meteo %>% 
  select(tmax_mean,tmin_mean, t5cm_min, t2m_mean_mon)
  • Macierz wykresów rozrzutu
pairs(~., data=sel, main="Macierz wykresów rozrzutu")

Macierz wykresów rozrzutu można także wykonać wykorzystując funkcję chart.Correlation() z pakietu PerformanceAnalytics

library("PerformanceAnalytics")
chart.Correlation(sel, histogram=TRUE, pch=19, method = "spearman" )
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties

  • Obliczenie macierzy współczynników korelacji dla wielu zmiennych.
cor(sel, use = "complete.obs", method = "spearman")
             tmax_mean tmin_mean  t5cm_min t2m_mean_mon
tmax_mean    1.0000000 0.8890966 0.7267208    0.9559309
tmin_mean    0.8890966 1.0000000 0.8723577    0.9395978
t5cm_min     0.7267208 0.8723577 1.0000000    0.8056521
t2m_mean_mon 0.9559309 0.9395978 0.8056521    1.0000000

Macierz korelacji wraz z informacją o poziomie istotności można uzyskać wykorzystując funckę rcorr() z pakietu Hmisc. 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). Wartość 0 oznacza, że p < 0.0001.

library(Hmisc)
kor <- rcorr(as.matrix(sel), type = "spearman")
print( kor, digits=3)
             tmax_mean tmin_mean t5cm_min t2m_mean_mon
tmax_mean         1.00      0.89     0.73         0.96
tmin_mean         0.89      1.00     0.87         0.94
t5cm_min          0.73      0.87     1.00         0.81
t2m_mean_mon      0.96      0.94     0.81         1.00

n= 708 


P
             tmax_mean tmin_mean t5cm_min t2m_mean_mon
tmax_mean               0         0        0          
tmin_mean     0                   0        0          
t5cm_min      0         0                  0          
t2m_mean_mon  0         0         0

6 Wizualizacja korelacji między zmiennymi z wykorzystaniem pakietu corrplot

Funkcja corrplot() z pakietu corrplot pozwala na wizualizację macierzy korelacji. Podstawowym parametrem jest macierz korelacji, którą można obliczyć wykorzystując funkcję cor() lub rcorr() z pakietu Hmisc.

library("corrplot")
corrplot 0.95 loaded
korelacja = cor(sel)
corrplot(korelacja)

library("corrplot")
library(Hmisc)
kor <- rcorr(as.matrix(sel), type = "spearman")
corrplot(kor$r)

library("corrplot")
library(Hmisc)
kor <- rcorr(as.matrix(sel), type = "spearman")
corrplot(kor$r, method = "number")

corrplot(kor$r, method = 'color', order = 'alphabet')

corrplot(kor$r, method = 'color', order = 'hclust')

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

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

corrplot.mixed(kor$r, lower = 'shade', upper = 'pie', order = 'hclust')

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

7 Wizualizacja macierzy wykresów rozrzutu

  • Funckja scatterplotMatrix() z pakietu car
library(car)
scatterplotMatrix(~tmax_mean+tmin_mean+t2m_mean_mon+t5cm_min, data=sel)

  • Funkcja chart.Correlation() z pakietu PerformanceAnalytics
library(PerformanceAnalytics)
chart.Correlation(sel, histogram=TRUE, pch=19, method = "spearman")
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties
Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
Cannot compute exact p-value with ties

  • Funkcja pairs.panels() z pakietu psych
library(psych)
pairs.panels(sel, scale=TRUE, method="spearman")