Výpočet korelace v R

Autor: Lenka Fiřtová

V tomto článku je popsán způsob, jak v R spočítat korelační koeficient mezi dvěma veličinami, respektive korelační matici v případě více veličin. Pro připomenutí, co to vlastně korelace je, viz článek s vysvětlením.

Výpočet korelačního koeficientu

Pro výpočty korelace použijeme  dataset trees, který je přímo obsažen v R. Tento dataset obsahuje pro 31 stromů údaje o jejich průměru (proměnná girth), výšce (proměnná height) a objemu dřeva (proměnná volume). Podívejme se, jak vypadá prvních pár řádků datasetu:

> head(trees)
   Girth  Height  Volume
1    8.3    70     10.3 
2    8.6    65     10.3
3    8.8    63     10.2
4   10.5    72     16.4
5   10.7    81     18.8
6   10.8    83     19.7

K výpočtu korelace slouží funkce cor.

Chceme-li vypočítat korelaci mezi dvěma veličinami (například průměrem a výškou stromů), zadáme funkci tyto dvě veličiny (syntaxe je: název datasetu $ název veličiny, případně název datasetu[ , číslo sloupce]). Nic dalšího není potřeba. Funkce spočítá tzv. Pearsonův korelační koeficient, což je právě ten korelační koeficient, který máme obvykle na mysli, když mluvíme o „korelaci“. Jedná se o podíl kovariance veličin a součinu jejich směrodatných odchylek (funkce dále umí spočítat Spearmanův koeficient korelace pořadí a Kendallův korelační koeficient, ty však nejsou předmětem tohoto článku).

> cor(trees$Girth, trees$Height)

[1] 0.5192801

Případně:

> cor(trees[ , 1], trees[ ,2])

[1] 0.5192801

Jak bychom očekávali, korelace je kladná – stromy s větším průměrem jsou i vyšší.

Problém může nastat ve chvíli, kdy data obsahují chybějící hodnoty (NA). Vytvoříme si dataset trees2, kde přidáme příkazem rbind k původnímu datasetu trees další řádek s chybějící informací o průměru (zbylé hodnoty si náhodně vymyslíme).

> trees2 = rbind(trees, c(NA, 90, 60))

Funkce cor pak vrátí chybu.

> cor(trees2$Girth, trees2$Height)

[1] NA

V případě chybějících pozorování proto musíme funkci říct, že má pracovat pouze se stromy, kde jsou všechny údaje k dispozici. K tomu slouží argument use = „complete.obs“.

> cor(trees2$Girth, trees2$Height, use = "complete.obs")

[1] 0.5192801

 

Korelační matice

Vraťme se nyní k původnímu datasetu. Rádi bychom zobrazili korelaci každé veličiny s každou najednou. V takovém případě stačí funkci cor dát více číselných veličin naráz, v našem případě klidně celý dataset trees, protože žádné nečíselné veličiny neobsahuje.

> cor(trees)
 
         Girth    Height    Volume
Girth  1.0000000 0.5192801 0.9671194 
Height 0.5192801 1.0000000 0.5982497
Volume 0.9671194 0.5982497 1.0000000

Kdyby dataset obsahoval například ještě čtvrtý sloupec s umístěním stromu (nečíselnou veličinou), řekli bychom funkci cor, ať počítá pouze s prvními třemi sloupci, následovně:

> cor(trees[,1:3])

Funkce cor v případě více než dvou veličin vrátí tzv. korelační matici. Na její diagonále jsou jedničky (tzn. jak koreluje veličina sama se sebou), mimo diagonálu pak korelace dílčích dvojic. Matice je symetrická (nad i pod diagonálou najdeme na odpovídajících pozicích stejná čísla).

Z korelační matice například na první pohled vidíme, že objem dřeva (volume) souvisí mnohem více s průměrem stromu (girth, korelace 0,97) než s výškou stromu (height, korelace 0,598). 

 

Test významnosti korelace

Chceme-li otestovat, jestli je korelace mezi veličinami statisticky významná, použijeme funkci cor.test.  Testem významnosti korelace zkoumáme, jestli korelace jen „zrovna nějak vyšla“ v našem konkrétním vzorku stromů, nebo jestli můžeme učinit závěr, že skutečně existuje nenulová korelace mezi veličinami v celé populaci. Podrobněji viz obecný článek o korelaci.

Prozkoumáme významnost korelace mezi proměnnými Girth a Height.

> cor.test(trees$Girth, trees$Height)

R vrátí následující výstup. Pod t se skrývá spočítaná testová statistika (3,2722) – tu bychom mohli porovnat s kritickou hodnotou z příslušných tabulek. Nemusíme ji však hledat, R totiž zobrazuje i tzv. p-value, česky p-hodnotu. K učinění závěru stačí porovnat p-hodnotu se zvolenou hladinou významnosti. Hladina významnosti činí obvykle 0,05. Pokud je p-hodnota menší než 0,05, což v našem případě je, pak můžeme tvrdit, že mezi veličinami skutečně existuje statisticky významný vztah.

data: trees$Girth and trees$Height
t = 3.2722, df = 29, p-value = 0.002758
alternative hypothesis: true correlation is not equal to 0
95 confidence interval:
0.2021327   0.7378538
sample estimates
cor
0.5192801

Vidíme, že 0.002758 je menší než 0.05. Mezi průměrem a výškou stromů existuje statisticky významná korelace.