Macierz korelacji Pearsona i Spearmana w R z przykładem

Korelacja dwuwymiarowa w R

Relacja dwuwymiarowa opisuje relację lub korelację pomiędzy dwiema zmiennymi w R. W tym samouczku omówimy koncepcję korelacji i pokażemy, jak można ją wykorzystać do pomiaru związku pomiędzy dowolnymi dwiema zmiennymi w R.

Korelacja w programowaniu R

Istnieją dwie podstawowe metody obliczania korelacji między dwiema zmiennymi w programowaniu w języku R:

• Pearson: Korelacja parametryczna
• Włócznik: Korelacja nieparametryczna

Macierz korelacji Pearsona w R

Metoda korelacji Pearsona jest zwykle stosowana jako podstawowe sprawdzenie związku między dwiema zmiennymi.

Współczynnik korelacji, napisane r, mierzy siłę liniowy związek między dwiema zmiennymi x i y. Oblicza się go w następujący sposób:

w

• jest odchyleniem standardowym x
• jest odchyleniem standardowym y

Korelacja mieści się w przedziale od -1 do 1.

• Wartość r bliska lub równa 0 oznacza niewielką lub żadną zależność liniową pomiędzy x i y.
• Im r jest bliższe 1 lub -1, tym silniejsza jest zależność liniowa.

Możesz sprawdzić, czy r różni się od zera za pomocą poniższej statystyki t, porównując ją z rozkładem Studenta z n – 2 stopniami swobody:

Korelacja rangi Spearmana w R

Korelacja rangowa sortuje obserwacje według rangi i oblicza stopień podobieństwa między rangami. Korelacja rangowa ma tę zaletę, że jest odporna na obserwacje odstające i nie jest powiązana z rozkładem danych. Korelacja rangowa jest również właściwym wyborem dla zmiennych porządkowych.

Korelacja rang Spearmana, zapisywana jako rho, również mieści się w przedziale od -1 do 1, a wartości bliskie skrajności wskazują na silną zależność monotoniczną. Oblicza się ją w następujący sposób:

Licznik jest kowariancją między rangami x i y, a mianownik jest iloczynem ich odchyleń standardowych.

W R oba te parametry oblicza się za pomocą funkcji cor(), która przyjmuje trzy argumenty: x, y i metodę.

cor(x, y, method)

Argumenty:

• x: Pierwszy wektor
• y: Drugi wektor
• metoda: Wzór używany do obliczenia korelacji. Trzy wartości ciągów:
  • "osoba"
  • „kendall”
  • "włócznik"

Można dodać opcjonalny argument, jeśli wektory zawierają brakującą wartość: use = „complete.obs”

Będziemy korzystać ze zbioru danych BudgetUK. Ten zbiór danych przedstawia alokację budżetu brytyjskich gospodarstw domowych w latach 1980–1982. Zawiera 1519 obserwacji obejmujących dziesięć cech, a wśród nich:

• jedzenie: podziel się jedzeniem, podziel się wydatkami
• wpaliwo: podziel się wydatkami na paliwo
• tkanina: udział w budżecie na wydatki na odzież
• walc: podziel się wydatkami na alkohol
• wtrans: podziel się wydatkami na transport
• kobieta: udział w wydatkach na inne towary
• toteks: całkowite wydatki gospodarstwa domowego w funtach
• dochód: całkowity dochód netto gospodarstwa domowego
• wiek: wiek gospodarstwa domowego
• dzieci: Liczba dzieci

Przykład

library(dplyr)
PATH <- "https://raw.githubusercontent.com/guru99-edu/R-Programming/master/british_household.csv"
data <- read.csv(PATH) %>%
    filter(income < 500) %>%
    mutate(log_income = log(income),
           log_totexp = log(totexp),
           children_fac = factor(children, order = TRUE, labels = c("No", "Yes"))) %>%
    select(-c(X, X.1, children, totexp, income))
glimpse(data)

Code Wyjaśnienie

• Najpierw importujemy dane i przyglądamy się funkcji przebłysk() z biblioteki dplyr.
• Trzy gospodarstwa domowe zgłosiły dochód wynoszący 500 lub więcej, więc filtr(dochód < 500) usuwa je, a liczba wierszy spada z 1,519 do 1,516.
• Powszechną praktyką jest przeliczanie zmiennej pieniężnej w logu. Pomaga zmniejszyć wpływ wartości odstających i zmniejsza skośność zbioru danych.

Wyjście:

## Observations: 1,516
## Variables: 10
## $ wfood        <dbl> 0.4272, 0.3739, 0.1941, 0.4438, 0.3331, 0.3752, 0...
## $ wfuel        <dbl> 0.1342, 0.1686, 0.4056, 0.1258, 0.0824, 0.0481, 0...
## $ wcloth       <dbl> 0.0000, 0.0091, 0.0012, 0.0539, 0.0399, 0.1170, 0...
## $ walc         <dbl> 0.0106, 0.0825, 0.0513, 0.0397, 0.1571, 0.0210, 0...
## $ wtrans       <dbl> 0.1458, 0.1215, 0.2063, 0.0652, 0.2403, 0.0955, 0...
## $ wother       <dbl> 0.2822, 0.2444, 0.1415, 0.2716, 0.1473, 0.3431, 0...
## $ age          <int> 25, 39, 47, 33, 31, 24, 46, 25, 30, 41, 48, 24, 2...
## $ log_income   <dbl> 4.867534, 5.010635, 5.438079, 4.605170, 4.605170,...
## $ log_totexp   <dbl> 3.912023, 4.499810, 5.192957, 4.382027, 4.499810,...
## $ children_fac <ord> Yes, Yes, Yes, Yes, No, No, No, No, No, No, Yes, ...

Współczynnik korelacji pomiędzy zmiennymi dochodu i wfood możemy obliczyć metodami „pearsona” i „spearmana”.

cor(data$log_income, data$wfood, method = "pearson")

Wyjście:

## [1] -0.2466986

cor(data$log_income, data$wfood, method = "spearman")

Wyjście:

## [1] -0.2501252

Zanim rozszerzymy to na każdą parę zmiennych, warto ustalić, jak należy odczytywać pojedynczy współczynnik.

Jak interpretować współczynnik korelacji

Współczynnik jest przydatny tylko wtedy, gdy można go określić. Poniższe paski to odczyt konwencjonalny, a znak jest odczytywany niezależnie od siły.

Wartość -0.2467 obliczona wcześniej pomiędzy log_income i wfood wskazuje zatem na słabą ujemną zależność: bogatsze gospodarstwa domowe przeznaczają nieco mniejszą część swojego budżetu na żywność.

Do każdego współczynnika odnoszą się trzy ostrzeżenia.

• Korelacja nie jest przyczynowością. Silne r oznacza, że ​​dwie zmienne poruszają się razem, a nie, że jedna powoduje drugą. Trzecia, niezmierzona zmienna często napędza obie.
• Pearson widzi tylko linie proste. Idealna relacja w kształcie litery U zwraca r bliskie zeru. Zawsze przed zaufaniem do liczby należy nanieść dane na wykres.
• Rozmiar ma większe znaczenie niż znaczenie. Przy 1,516 obserwacjach współczynnik 0.06 może okazać się statystycznie istotny, choć w praktyce bez znaczenia.

Jak testować istotność korelacji za pomocą cor.test()

Funkcja cor() zwraca współczynnik i nic więcej. Dla pojedynczej pary funkcja cor.test() dodaje wartość p i przedział ufności w jednym wywołaniu.

cor.test(data$log_income, data$wfood, method = "pearson")

Wynik składa się z czterech części wartych przeczytania.

1. t i df: statystyka testowa i jej stopnie swobody, n – 2.
2. Wartość p:prawdopodobieństwo uzyskania tak dużego współczynnika, jeśli prawdziwa korelacja wynosi zero.
3. 95-procentowy przedział ufności: prawdopodobny zakres prawdziwej korelacji. Jeśli wyklucza zero, zależność jest istotna na tym poziomie.
4. przykładowy szacunek: współczynnik sam w sobie, identyczny z tym, który zwraca funkcja cor().

Ta sama funkcja uruchamia testy oparte na rangach poprzez zmianę jednego argumentu:

# Spearman rank correlation with a p-value
cor.test(data$log_income, data$wfood, method = "spearman")

# One-sided test: is the correlation greater than zero?
cor.test(data$log_income, data$wfood, alternative = "greater")

Kiedy używać którego. Użyj funkcji cor.test(), gdy badasz jedną konkretną parę, ponieważ daje ona przedział ufności pomijany przez rcorr(). Użyj rcorr() z Hmisc, pokazanego powyżej, gdy potrzebujesz wartości p dla całej macierzy naraz. Należy pamiętać, że testowanie wielu par zawyża wskaźnik wyników fałszywie dodatnich, dlatego przed wyciągnięciem wniosków z dużej macierzy należy skorygować wartości p za pomocą p.adjust(p_value, method = “BH”).

Macierz korelacji w R

Korelacja dwuwymiarowa to dobry początek, ale spojrzenie wielowymiarowe daje szerszy obraz. macierz korelacji jest kwadratową tabelą przechowującą parami korelację każdej zmiennej względem każdej innej.

Funkcja cor() zwraca macierz korelacji. Jedyna różnica w przypadku korelacji dwuwymiarowej polega na tym, że nie musimy określać, które zmienne. Domyślnie R oblicza korelację pomiędzy wszystkimi zmiennymi.

Nie można obliczyć korelacji dla czynnika, dlatego usuń każdą kolumnę kategoryjną przed przekazaniem ramki danych do cor().

Macierz korelacji jest symetryczna, co oznacza, że ​​wartości powyżej przekątnej mają takie same wartości jak ta poniżej. Bardziej wizualne jest pokazanie połowy macierzy.

children_fac jest wykluczone, ponieważ cor() nie może działać na czynniku.

# the last column of data is a factor level. We don't include it in the code
mat_1 <-as.dist(round(cor(data[,1:9]),2))
mat_1

Code Wyjaśnienie

• cor(dane[, 1:9]):Oblicz macierz korelacji na podstawie dziewięciu kolumn liczbowych
• okrągły(…, 2):Zaokrąglij każdy współczynnik do dwóch miejsc po przecinku
• as.dist():Wydrukuj tylko dolny trójkąt, ponieważ macierz jest symetryczna

Wyjście:

##            wfood wfuel wcloth  walc wtrans wother   age log_income
## wfuel       0.11                                                  
## wcloth     -0.33 -0.25                                            
## walc       -0.12 -0.13  -0.09                                     
## wtrans     -0.34 -0.16  -0.19 -0.22                               
## wother     -0.35 -0.14  -0.22 -0.12  -0.29                        
## age         0.02 -0.05   0.04 -0.14   0.03   0.02                 
## log_income -0.25 -0.12   0.10  0.04   0.06   0.13  0.23           
## log_totexp -0.50 -0.36   0.34  0.12   0.15   0.15  0.21       0.49

Poziom istotności

Sam współczynnik nie mówi, czy relacja jest statystycznie wiarygodna. Funkcja rcorr() z biblioteki Hmisc zwraca wartość p dla każdej pary. Bibliotekę można pobrać ze strony Conda i skopiuj kod, aby wkleić go w terminalu:

conda install -c r r-hmisc

Funkcja rcorr() wymaga przechowywania ramki danych w postaci macierzy. Możemy przekształcić nasze dane w macierz przed obliczeniem macierzy korelacji z wartością p.

library("Hmisc")
data_rcorr <-as.matrix(data[, 1: 9])

mat_2 <-rcorr(data_rcorr)
# mat_2 <-rcorr(as.matrix(data)) returns the same output

Obiekt listy mat_2 zawiera trzy elementy:

• r: Dane wyjściowe macierzy korelacji
• n: Liczba obserwacji
• P: wartość p

Nas interesuje trzeci element, wartość p. Często zamiast współczynnika korelacji przedstawia się macierz korelacji z wartością p.

p_value <-round(mat_2[["P"]], 3)
p_value

Code Wyjaśnienie

• mat_2[[„P”]]: Wartości p są przechowywane w elemencie zwanym P
• round(mat_2[[„P”]], 3): Zaokrąglij elementy trzema cyframi

Wyjście:

           wfood wfuel wcloth  walc wtrans wother   age log_income log_totexp
wfood         NA 0.000  0.000 0.000  0.000  0.000 0.365      0.000          0
wfuel      0.000    NA  0.000 0.000  0.000  0.000 0.076      0.000          0
wcloth     0.000 0.000     NA 0.001  0.000  0.000 0.160      0.000          0
walc       0.000 0.000  0.001    NA  0.000  0.000 0.000      0.105          0
wtrans     0.000 0.000  0.000 0.000     NA  0.000 0.259      0.020          0
wother     0.000 0.000  0.000 0.000  0.000     NA 0.355      0.000          0
age        0.365 0.076  0.160 0.000  0.259  0.355    NA      0.000          0
log_income 0.000 0.000  0.000 0.105  0.020  0.000 0.000         NA          0
log_totexp 0.000 0.000  0.000 0.000  0.000  0.000 0.000      0.000         NA

Wizualizacja macierzy korelacji w R

Mapa cieplna to kolejny sposób odczytu macierzy korelacji. Biblioteka GGally rozszerza ggplot2 i jest instalowana z CRAN, a nie z Conda:

install.packages("GGally")

Biblioteka zawiera różne funkcje umożliwiające pokazanie statystyk podsumowujących, takich jak korelacja i rozkład wszystkich zmiennych w a matryca.

Funkcja ggcorr() ma wiele argumentów. Przedstawimy tylko argumenty, których będziemy używać w tutorialu:

Funkcja ggcorr

ggcorr(df, method = c("pairwise", "pearson"),
  nbreaks = NULL, digits = 2, low = "#3B9AB2",
  mid = "#EEEEEE", high = "#F21A00",
  geom = "tile", label = FALSE,
  label_alpha = FALSE)

argumenty:

• df: Wykorzystano zbiór danych
• metoda: Wzór do obliczania korelacji. Domyślnie obliczane są parami i Pearson
• nprzerw: Zwraca kategoryczny zakres zabarwienia współczynników. Domyślnie nie ma przerwy, a gradient kolorów jest ciągły
• cyfry: Zaokrąglij współczynnik korelacji. Domyślnie ustaw na 2
• Niska: Sterowanie dolnym poziomem zabarwienia
• średni: Kontroluj środkowy poziom zabarwienia
• wysoka: Kontroluj wysoki poziom zabarwienia
• geom: Kontroluj kształt argumentu geometrycznego. Domyślnie „kafelek”
• etykieta: Wartość logiczna. Wyświetlaj etykietę lub nie. Domyślnie ustawione na „FAŁSZ”.

Podstawowa mapa cieplna

Najbardziej podstawowym wykresem pakietu jest mapa cieplna. Legenda wykresu pokazuje kolor gradientowy od – 1 do 1, przy czym kolor gorący oznacza silną korelację dodatnią, a kolor zimny korelację ujemną.

library(GGally)
ggcorr(data)

Code Wyjaśnienie

• ggcorr(dane): Potrzebny jest tylko jeden argument, którym jest nazwa ramki danych. Zmienne na poziomie współczynnika nie są uwzględnione na wykresie.

Wyjście:

Dodawanie kontroli do mapy cieplnej

Do wykresu możemy dodać więcej kontrolek:

ggcorr(data,
    nbreaks = 6,
    low = "steelblue",
    mid = "white",
    high = "darkred",
    geom = "circle")

Code Wyjaśnienie

• nprzerw=6: przełam legendę za pomocą 6 stopni.
• niski = „stalowoniebieski”: Użyj jaśniejszych kolorów, aby uzyskać ujemną korelację
• środek = „biały”: Użyj białych kolorów dla korelacji średnich zakresów
• wysoki = „ciemnoczerwony”: Użyj ciemnych kolorów, aby uzyskać dodatnią korelację
• geom = „okrąg”: Użyj okręgu jako kształtu okien na mapie cieplnej. Rozmiar okręgu jest proporcjonalny do wartości bezwzględnej korelacji.

Wyjście:

Dodawanie etykiety do mapy cieplnej

GGally pozwala nam dodać etykietę wewnątrz okien:

ggcorr(data,
    nbreaks = 6,
    label = TRUE,
    label_size = 3,
    color = "grey50")

Code Wyjaśnienie

• etykieta = PRAWDA: Dodaj wartości współczynników korelacji wewnątrz mapy cieplnej.
• kolor = „szary50”: Wybierz kolor, np. szary
• rozmiar_etykiety = 3: Ustaw rozmiar etykiety na 3

Wyjście:

Funkcja ggpairs

Biblioteka GGally udostępnia również funkcję ggpairs(), która zwraca macierz wykresów. Dla k wybranych zmiennych wynik to siatka k na k: przekątna pokazuje rozkład każdej zmiennej, a panele powyżej i poniżej przekątnej mogą przeprowadzać różne obliczenia. Składnia jest następująca:

ggpairs(df, columns = 1:ncol(df), title = NULL,
    upper = list(continuous = "cor"),
    lower = list(continuous = "smooth"),
    mapping = NULL)		

argumenty:

• df: Wykorzystano zbiór danych
• kolumny: Wybierz kolumny, aby narysować wykres
• tytuł: Podaj tytuł
• górny: Steruj polami nad przekątną wykresu. Należy podać typ obliczeń lub wykresu do zwrócenia. Jeśli ciągła = „cor”, prosimy R o obliczenie korelacji. Należy pamiętać, że argument musi być listą. Dostępne są inne argumenty; patrz Dokumentacja GGally po więcej informacji.
• niższy:Steruj polami znajdującymi się poniżej przekątnej.
• mapaping: Wskazuje estetykę wykresu. Na przykład możemy obliczyć wykres dla różnych grup.

Analiza dwuwymiarowa z ggpair z grupąping

Następny wykres przedstawia trzy informacje:

• Macierz korelacji między log_totex, log_income, wiekiem i zmienną wtrans pogrupowaną według tego, czy gospodarstwo domowe ma dziecko, czy nie.
• Narysuj rozkład każdej zmiennej według grup
• Wyświetl wykres punktowy z trendem według grup

library(ggplot2)
ggpairs(data, columns = c("log_totexp", "log_income", "age", "wtrans"), title = "Bivariate analysis of revenue expenditure by the British household", upper = list(continuous = wrap("cor",
        size = 3)),
    lower = list(continuous = wrap("smooth",
        alpha = 0.3,
        size = 0.1)),
    mapping = aes(color = children_fac))

Code Wyjaśnienie

• kolumny = c(„log_totex”, „log_income”, „wiek”, „wtrans”): wybierz zmienne, które mają być pokazane na wykresie
• title = „Dwuwymiarowa analiza wydatków dochodowych brytyjskiego gospodarstwa domowego”: Dodaj tytuł
• górna = lista(): Kontroluj górną część wykresu. Tj. Powyżej przekątnej
• ciągły = zawiń(„cor”, rozmiar = 3)): Oblicz współczynnik korelacji. Zawijamy argument Continuous wewnątrz funkcji wrap(), aby kontrolować estetykę wykresu (tj. rozmiar = 3). -lower = list(): Kontroluje dolną część wykresu. Tzn. poniżej przekątnej.
• ciągły = zawiń („gładki”, alfa = 0.3, rozmiar = 0.1): Dodaj wykres punktowy z trendem liniowym. Zawijamy argument ciągły wewnątrz funkcji wrap(), aby kontrolować estetykę wykresu (tj. rozmiar=0.1, alfa=0.3).
• mapaping = aes(kolor = dzieci_fac):Podziel każdy panel według children_fac, uporządkowanego czynnika oznaczonego jako „Nie” dla gospodarstw domowych bez dzieci i „Tak” dla gospodarstw domowych z dziećmi

Wyjście:

Analiza dwuwymiarowa z ggpair z grupami częściowymiping

Poniższy wykres jest nieco inny. Zmieniamy położenie mapy.ping wewnątrz górnego argumentu.

ggpairs(data, columns = c("log_totexp", "log_income", "age", "wtrans"),
    title = "Bivariate analysis of revenue expenditure by the British household",
    upper = list(continuous = wrap("cor",
            size = 3),
        mapping = aes(color = children_fac)),
    lower = list(
        continuous = wrap("smooth",
            alpha = 0.3,
            size = 0.1))
)

Code Wyjaśnienie

• Dokładny ten sam kod, co w poprzednim przykładzie, z wyjątkiem:
• mapaping = aes(color = children_fac): Przenieś listę w upper = list(). Chcemy, aby obliczenia były ułożone w stosy według grup tylko w górnej części grafu.

Wyjście:

Korelacja w R: najważniejsze wnioski i odniesienie do funkcji

• Relacja dwuwymiarowa opisuje związek – lub korelację – pomiędzy dwiema zmiennymi w R.
• Istnieją dwie podstawowe metody obliczania korelacji między dwiema zmiennymi w Programowanie R.: Pearson i Spearman.
• Metoda korelacji Pearsona jest zwykle stosowana jako podstawowe sprawdzenie związku między dwiema zmiennymi.
• Korelacja rang sortuje obserwacje według rang i oblicza poziom podobieństwa między rangami.
• Współczynnik korelacji rang Spearmana mieści się w przedziale od -1 do 1, a wartości bliskie obu skrajności wskazują na silną zależność monotoniczną.
• Macierz korelacji to kwadratowa tabela zawierająca parami korelację każdej zmiennej.
• Wartość p informuje, czy zaobserwowana korelacja jest statystycznie odróżnialna od zera.

Poniżej wymieniono wszystkie funkcje korelacji użyte w tym samouczku:

cor(dfx2, method = "pearson")

cor(dfx2, method = "spearman")

cor(df, method = "pearson")

cor(df, method = "spearman")

rcorr(as.matrix(data[,1:9]))[["P"]]

ggcorr(df)

ggpairs(df, columns = c("x1", "x2"))

FAQ