Pearsonova i Spearmanova korelacijska matrica u R-u s primjerom

Bivarijatna korelacija u R

Bivarijantni odnos opisuje odnos - ili korelaciju - između dviju varijabli u R. U ovom vodiču ćemo raspravljati o konceptu korelacije i pokazati kako se može koristiti za mjerenje odnosa između bilo koje dvije varijable u R.

Korelacija u R programiranju

Postoje dvije primarne metode za izračunavanje korelacije između dvije varijable u R programiranju:

• Pearson: Parametarska korelacija
• Kopljanik: Neparametarska korelacija

Pearsonova korelacijska matrica u R

Pearsonova metoda korelacije obično se koristi kao primarna provjera odnosa između dviju varijabli.

The koeficijent korelacije, označeno kao r, mjeri snagu linearan odnos između dvije varijable x i y. Izračunava se na sljedeći način:

sa

• je standardna devijacija x
• je standardna devijacija y

Korelacija se kreće između -1 i 1.

• Vrijednost r blizu ili jednaka 0 implicira malu ili nikakvu linearnu vezu između x i y.
• Što je r bliže 1 ili -1, to je linearni odnos jači.

Možete provjeriti razlikuje li se r od nule pomoću t-statistike u nastavku, uspoređujući je sa Studentovom distribucijom s n – 2 stupnja slobode:

Spearmanova korelacija ranga u R

Rang korelacija sortira opažanja po rangu i izračunava razinu sličnosti između rangova. Rang korelacija ima prednost što je robusna na outliere i nije povezana s distribucijom podataka. Rang korelacija je također pravi izbor za ordinalne varijable.

Spearmanova korelacija ranga, označena kao rho, također se kreće od -1 do 1, a vrijednosti blizu bilo kojem ekstremu ukazuju na jaku monotonu vezu. Izračunava se na sljedeći način:

Brojnik je kovarijanca između rangova x i y, a nazivnik je umnožak njihovih standardnih devijacija.

U R-u se oboje izračunava funkcijom cor(), koja uzima tri argumenta: x, y i metodu.

cor(x, y, method)

argumenti:

• x: Prvi vektor
• y: Drugi vektor
• način: Formula koja se koristi za izračunavanje korelacije. Tri vrijednosti niza:
  • “pearson”
  • “kendall”
  • "kopljanik"

Neobavezni argument može se dodati ako vektori sadrže vrijednost koja nedostaje: use = “complete.obs”

Koristit ćemo skup podataka BudgetUK. Ovaj skup podataka izvješćuje o raspodjeli proračuna britanskih kućanstava između 1980. i 1982. Postoji 1519 promatranja s deset značajki, među kojima su:

• wfood: dijeliti hranu dijeliti trošiti
• wgorivo: podijelite potrošnju goriva
• krpa: udio proračuna za potrošnju na odjeću
• hodati: podijeliti trošenje alkohola
• wtrans: podijelite troškove prijevoza
• majka: udio potrošnje ostalih dobara
• totexp: ukupna potrošnja kućanstva u funtama
• dohodak: ukupni neto prihod kućanstva
• starost: starost kućanstva
• djeca: broj djece

Primjer

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 Objašnjenje

• Prvo uvozimo podatke i pogledamo pomoću funkcije glimpse() iz knjižnice dplyr.
• Tri kućanstva prijavljuju prihod od 500 ili više, pa ih filter(prihod < 500) uklanja i broj redaka pada s 1,519 na 1,516.
• Uobičajena je praksa pretvaranje monetarne varijable u log. Pomaže smanjiti utjecaj odstupanja i smanjuje asimetriju u skupu podataka.

Izlaz:

## 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, ...

Možemo izračunati koeficijent korelacije između varijabli dohotka i wfood metodama “pearson” i “spearman”.

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

Izlaz:

## [1] -0.2466986

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

Izlaz:

## [1] -0.2501252

Prije proširenja ovoga na svaki par varijabli, vrijedi utvrditi kako treba čitati pojedinačni koeficijent.

Kako interpretirati koeficijent korelacije

Koeficijent je koristan samo kada možete reći što znači. Trake ispod su konvencionalno očitanje, a predznak se očitava odvojeno od jačine.

Vrijednost od -0.2467 izračunata ranije između log_income i wfood stoga je slaba negativna veza: bogatija kućanstva troše nešto manji dio svog proračuna na hranu.

Tri upozorenja primjenjuju se na svaki koeficijent.

• Korelacija nije uzročnost. Snažan r kaže da se dvije varijable kreću zajedno, nikada da jedna uzrokuje drugu. Treća, neizmjerena varijabla često pokreće obje.
• Pearson vidi samo ravne linije. Savršen odnos u obliku slova U vraća r blizu nule. Uvijek iscrtajte podatke prije nego što se povjerujete broju.
• Veličina je važnija od značaja. S 1,516 opažanja, koeficijent od 0.06 može biti statistički značajan, a ipak praktički beznačajan.

Kako testirati značajnost korelacije pomoću cor.test()

cor() vraća koeficijent i ništa drugo. Za jedan par, cor.test() zbraja p-vrijednost i interval pouzdanosti u jednom pozivu.

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

Izlaz ima četiri dijela vrijedna čitanja.

1. t i df: testna statistika i njezini stupnjevi slobode, n – 2.
2. p-vrijednostvjerojatnost da će koeficijent biti ovako velik ako bi stvarna korelacija bila nula.
3. 95-postotni interval pouzdanosti: vjerojatni raspon za stvarnu korelaciju. Ako isključuje nulu, odnos je značajan na toj razini.
4. primjer procjene: sam koeficijent, identičan onome što vraća cor().

Ista funkcija pokreće testove temeljene na rangu promjenom jednog argumenta:

# 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")

Kada koji koristiti. Koristite cor.test() kada ispitujete jedan određeni par, jer daje interval pouzdanosti koji rcorr() izostavlja. Koristite rcorr() iz Hmisc-a, prikazan gore, kada trebate p-vrijednosti za cijelu matricu odjednom. Imajte na umu da testiranje mnogih parova povećava stopu lažno pozitivnih rezultata, stoga prilagodite p-vrijednosti pomoću p.adjust(p_value, method = “BH”) prije izvlačenja zaključaka iz velike matrice.

Korelacijska matrica u R

Bivarijantna korelacija je dobar početak, ali multivarijantni pogled daje širu sliku. korelacijska matrica je kvadratna tablica koja sadrži parne korelacije svake varijable u odnosu na sve ostale.

Funkcija cor() vraća korelacijsku matricu. Jedina razlika u odnosu na bivarijatnu korelaciju je što ne trebamo specificirati koje varijable. Prema zadanim postavkama, R izračunava korelaciju između svih varijabli.

Korelacija se ne može izračunati za faktor, stoga izostavite svaki kategorički stupac prije nego što proslijedite podatkovni okvir funkciji cor().

Korelacijska matrica je simetrična što znači da vrijednosti iznad dijagonale imaju iste vrijednosti kao one ispod. Vizualnije je prikazati polovicu matrice.

children_fac je isključen jer cor() ne može operirati na faktoru.

# 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 Objašnjenje

• cor(podaci[, 1:9])Izračunajte matricu korelacije na devet numeričkih stupaca
• okruglo(…, 2)Zaokružite svaki koeficijent na dvije decimale
• as.dist()Ispišite samo donji trokut, budući da je matrica simetrična

Izlaz:

##            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

Razina značajnosti

Koeficijent sam po sebi ne govori je li odnos statistički pouzdan. Funkcija rcorr() iz Hmisc biblioteke vraća p-vrijednost za svaki par. Biblioteku možemo preuzeti s konda i kopirajte kod da biste ga zalijepili u terminal:

conda install -c r r-hmisc

Rcorr() zahtijeva da okvir podataka bude pohranjen kao matrica. Možemo pretvoriti naše podatke u matricu prije nego što izračunamo korelacijsku matricu s p-vrijednošću.

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

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

Objekt liste mat_2 sadrži tri elementa:

• r: Izlaz korelacijske matrice
• n: Broj opažanja
• P: p-vrijednost

Zanima nas treći element, p-vrijednost. Uobičajeno je prikazati korelacijsku matricu s p-vrijednošću umjesto koeficijenta korelacije.

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

Code Objašnjenje

• mat_2[[“P”]]: p-vrijednosti su pohranjene u elementu koji se zove P
• okruglo(mat_2[[“P”]], 3): Zaokružite elemente s tri znamenke

Izlaz:

           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

Vizualizacija korelacijske matrice u R

Toplinska karta je još jedan način čitanja korelacijske matrice. Biblioteka GGally proširuje ggplot2 i instalira se iz CRAN-a, a ne iz conda-e:

install.packages("GGally")

Knjižnica uključuje različite funkcije za prikaz sumarne statistike kao što je korelacija i distribucija svih varijabli u matrica.

Funkcija ggcorr() ima mnogo argumenata. Predstavit ćemo samo argumente koje ćemo koristiti u uputama:

Funkcija ggcorr

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

argumenti:

• df: Korišteni skup podataka
• način: Formula za izračunavanje korelacije. Prema zadanim postavkama izračunavaju se parovi i Pearson
• nbreaks: Vrati kategorički raspon za bojenje koeficijenata. Prema zadanim postavkama, nema prekida i gradijent boja je kontinuiran
• znamenki: Zaokružite koeficijent korelacije. Prema zadanim postavkama postavljeno na 2
• nisko: Kontrolirajte nižu razinu obojenosti
• srednji: Kontrolirajte srednju razinu obojenosti
• visok: Kontrolirajte visoku razinu obojenosti
• geom: Kontrolirajte oblik geometrijskog argumenta. Prema zadanim postavkama, "pločica"
• oznaka: Booleova vrijednost. Prikaži ili ne oznaku. Prema zadanim postavkama postavljeno na `FALSE`

Osnovna toplinska karta

Najosnovniji crtež paketa je toplinska karta. Legenda grafikona prikazuje gradijent boje od – 1 do 1, s vrućom bojom koja označava jaku pozitivnu korelaciju, a hladnom bojom negativnu korelaciju.

library(GGally)
ggcorr(data)

Code Objašnjenje

• ggcorr(podaci): Potreban je samo jedan argument, a to je naziv okvira podataka. Varijable razine faktora nisu uključene u dijagram.

Izlaz:

Dodavanje kontrole na toplinsku kartu

Grafikonu možemo dodati više kontrola:

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

Code Objašnjenje

• nprekidi=6: razbiti legendu sa 6 redova.
• niska = "čeličnoplava": Koristite svjetlije boje za negativnu korelaciju
• sredina = “bijelo”: Koristite bijele boje za korelaciju srednjih raspona
• visoko = "tamno": Koristite tamne boje za pozitivnu korelaciju
• geom = “krug”: Koristite krug kao oblik prozora na toplinskoj karti. Veličina kruga proporcionalna je apsolutnoj vrijednosti korelacije.

Izlaz:

Dodavanje oznake na toplinsku kartu

GGally nam omogućuje dodavanje oznake unutar prozora:

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

Code Objašnjenje

• oznaka = TRUE: Dodajte vrijednosti koeficijenata korelacije unutar toplinske karte.
• boja = “siva50”: Odaberite boju, tj. sivu
• veličina_oznake = 3: Postavite veličinu naljepnice na 3

Izlaz:

Funkcija ggpairs

Biblioteka GGally također nudi ggpairs(), koja vraća matricu grafova. Za k odabranih varijabli rezultat je mreža ak x k: dijagonala prikazuje distribuciju svake varijable, dok ploče iznad i ispod dijagonale mogu nositi različit izračun. Sintaksa je:

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

argumenti:

• df: Korišteni skup podataka
• stupovi: Odaberite stupce za crtanje dijagrama
• naslov: Uključite naslov
• gornji: Kontrolirajte okvire iznad dijagonale grafikona. Potrebno je navesti vrstu izračuna ili grafa koji treba vratiti. Ako je continuous = "cor", tražimo od R-a da izračuna korelaciju. Imajte na umu da argument mora biti popis. Dostupni su i drugi argumenti; pogledajte GGally dokumentacija za više informacija.
• smanjiti: Kontrolirajte okvire ispod dijagonale.
• kartaping: Označava estetiku grafikona. Na primjer, možemo izračunati grafikon za različite skupine.

Bivarijantna analiza s ggpairom i grouping

Sljedeći grafikon prikazuje tri informacije:

• Korelacijska matrica između varijabli log_totex, log_income, age i wtrans grupiranih prema tome ima li kućanstvo dijete ili ne.
• Nacrtajte distribuciju svake varijable po skupini
• Prikažite dijagram raspršenosti s trendom po grupi

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 Objašnjenje

• stupci = c(“log_totex”, “log_income”, “age”, “wtrans”): Odaberite varijable za prikaz na grafikonu
• title = “Bivarijantna analiza rashoda prihoda britanskog kućanstva”: Dodajte naslov
• gornji = lista(): Kontrolirajte gornji dio grafikona. Tj. Iznad dijagonale
• kontinuirano = omot ("cor", veličina = 3)): Izračunajte koeficijent korelacije. Umotavamo kontinuirani argument unutar funkcije wrap() kako bismo kontrolirali estetiku grafa (tj. veličina = 3) -lower = list(): Kontrolirajte donji dio grafa. Tj. Ispod dijagonale.
• kontinuirano = omot (“glatko”, alfa = 0.3, veličina = 0.1): Dodajte dijagram raspršenosti s linearnim trendom. Umotavamo kontinuirani argument unutar funkcije wrap() kako bismo kontrolirali estetiku grafikona (tj. veličina=0.1, alfa=0.3)
• kartaping = aes(boja = dječji_fac)Podijelite svaki panel prema children_fac, uređenom faktoru označenom s "Ne" za kućanstva bez djece i "Da" za kućanstva s djecom

Izlaz:

Bivarijatna analiza s ggpairom i parcijalnom grupomping

Grafikon ispod je malo drugačiji. Mijenjamo položaj karteping unutar gornjeg argumenta.

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 Objašnjenje

• Potpuno isti kod kao prethodni primjer osim:
• kartaping = aes(color = children_fac): Premjesti popis u gornji = list(). Želimo da izračun bude složen po grupama samo u gornjem dijelu grafa.

Izlaz:

Korelacija u R-u: Ključne informacije i referenca funkcija

• Bivarijantni odnos opisuje odnos - ili korelaciju - između dvije varijable u R.
• Postoje dvije primarne metode za izračunavanje korelacije između dvije varijable R Programiranje: Pearson & Spearman.
• Pearsonova metoda korelacije obično se koristi kao primarna provjera odnosa između dviju varijabli.
• Korelacija ranga razvrstava zapažanja prema rangu i izračunava razinu sličnosti između ranga.
• Spearmanova korelacija ranga kreće se od -1 do 1, a vrijednosti blizu bilo kojeg ekstrema ukazuju na jaku monotonu vezu.
• Matrica korelacije je kvadratna tablica koja sadrži parne korelacije svake varijable.
• P-vrijednost vam govori je li uočena korelacija statistički različita od nule.

Svaka korelacijska funkcija korištena u ovom tutorijalu navedena je u nastavku:

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

Pitanja i odgovori