Korelacja w R: cor(), cor.test() i macierz korelacji

Korelacja między dwiema zmiennymi: cor()

Korelacja pyta: gdy jedna zmienna rośnie, czy druga też ma tendencję do wzrostu (dodatnia), do spadku (ujemna), czy robi swoje (blisko zera)? W R to jedno wywołanie:

print(cor(mtcars$wt, mtcars$mpg))

Wynik to około −0.87: cięższe samochody mają gorsze spalanie, a zależność jest silna. Ta jedna liczba to współczynnik korelacji Pearsona, powszechnie oznaczany r.

Jak czytać r

Współczynnik zawsze leży między −1 a +1. Znak podaje kierunek, wartość bezwzględna siłę:

| |r| | Typowa interpretacja | | --- | --- | | 0.0 do 0.2 | pomijalna | | 0.2 do 0.4 | słaba | | 0.4 do 0.6 | umiarkowana | | 0.6 do 0.8 | silna | | 0.8 do 1.0 | bardzo silna |

Traktuj te przedziały jako punkt wyjścia do rozmowy, a nie prawo: w fizyce r równe 0.6 rozczarowuje, w psychologii to szczyt kariery. Dwie właściwości warto przyswoić: r nie ma jednostek (korelacja masy w tonach z mpg daje to samo r co masa w kilogramach z mpg, bo r liczy się na wartościach standaryzowanych), a r mierzy tylko zależność liniową, więc idealna zależność w kształcie litery U może mieć r ≈ 0.

I zdanie, które trzeba powiedzieć: korelacja to nie przyczynowość. Sprzedaż lodów i liczba utonięć silnie korelują w kolejnych miesiącach roku, ale nie dlatego, że lody topią ludzi, tylko dlatego, że lato napędza jedno i drugie. Korelacja mówi, że dwie zmienne zmieniają się razem; milczy o tym, dlaczego. Może x wpływa na y, może y na x, a może trzecia rzecz (tutaj pora roku) wpływa na obie. Rozstrzygnięcie tego wymaga eksperymentów albo starannego rozumowania przyczynowego, a nie większego r.

Spearman i Kendall: kiedy Pearson to złe narzędzie

Pearson działa na surowych wartościach, więc jest wrażliwy na wartości odstające i ślepy na zależności zakrzywione. Argument method przełącza na alternatywy oparte na rangach:

x <- c(1, 2, 3, 4, 5, 6, 7, 8)
y <- c(1, 4, 9, 16, 25, 36, 49, 300)  # increasing, with a wild last value

print(cor(x, y))                        # Pearson - dragged by the outlier
print(cor(x, y, method = "spearman"))   # ranks: perfectly monotonic
print(cor(x, y, method = "kendall"))

Spearman zastępuje każdą wartość jej rangą, a potem liczy Pearsona na rangach. Ponieważ y zawsze rośnie, gdy rośnie x, wszystkie rangi się zgadzają i Spearman zwraca dokładnie 1: wielkość wartości odstającej przestaje mieć znaczenie, liczy się tylko jej pozycja. Sięgnij po Spearmana, gdy dane są porządkowe (skale ankietowe), mocno skośne albo gdy zależność jest monotoniczna, ale nie prosta. Kendall odpowiada na podobne pytanie na podstawie par zgodnych i niezgodnych; jest odporniejszy przy małych próbach, ale wolniejszy, a domyślnym wyborem zwykle jest Spearman.

Macierz korelacji

Żeby przejrzeć zależności między wieloma zmiennymi naraz, przekaż do cor() kilka kolumn liczbowych:

m <- cor(mtcars[, c("mpg", "wt", "hp", "disp")])
print(round(m, 2))

Każda zmienna względem każdej innej, z jedynkami na przekątnej (wszystko koreluje idealnie samo ze sobą) i lustrzanym odbiciem po obu jej stronach. Zaokrąglenie do dwóch miejsc po przecinku ma większe znaczenie, niż się wydaje: niezaokrąglona macierz to ściana cyfr, a macierz ma się dać przejrzeć. Tutaj widać, że mpg koreluje ujemnie ze wszystkimi trzema (cięższe, mocniejsze auta z większym silnikiem palą więcej), a wt, hp i disp są ze sobą silnie dodatnio skorelowane. To grupa zmiennych "duże auto", która będzie ważna przy regresji liniowej i jej problemach ze współliniowością.

Brakujące wartości: argument use

Przy brakujących danych cor() domyślnie zwraca NA, zamiast zgadywać:

d <- mtcars[, c("mpg", "wt", "hp")]
d$wt[c(3, 7)] <- NA  # knock two values out

print(cor(d)[1, 2])
print(cor(d, use = "complete.obs")[1, 2])
print(cor(d, use = "pairwise.complete.obs")[1, 2])

• use = "complete.obs" usuwa każdy wiersz zawierający jakiekolwiek NA, a potem liczy całą macierz z pozostałych wierszy. Spójnie, ale rozrzutnie: brakujące wt usuwa ten wiersz także z pary mpg i hp.
• use = "pairwise.complete.obs" liczy każdą komórkę ze wszystkich wierszy, w których obecna jest ta konkretna para. Zachowuje więcej danych, ale różne komórki opierają się na różnych podzbiorach, co bardzo rzadko może dać macierz wewnętrznie niespójną.

Przy kilku zabłąkanych NA każda opcja jest w porządku; po prostu napisz, której użyto.

Czy to istotne? cor.test()

cor() daje liczbę, ale nie mówi, czy może to być szum. cor.test() przeprowadza test hipotezy:

print(cor.test(mtcars$wt, mtcars$mpg))

Przejdźmy przez wynik blok po bloku:

• t = −9.56, df = 30: statystyka testowa. Hipoteza zerowa mówi, że prawdziwa korelacja wynosi zero; obserwowane r jest przeliczane na statystykę t z n − 2 stopniami swobody (32 samochody − 2).
• p-value = 1.29e-10: gdyby prawdziwa korelacja wynosiła zero, prawdopodobieństwo zobaczenia r tak dalekiego od zera w próbie 32 aut to około 0.0000000001. To przytłaczający dowód, że zależność jest prawdziwa. Pamiętaj jednak, że p-value mówi o tym, czy r różni się od zera, a nie o tym, czy zależność jest duża albo przyczynowa.
• 95 percent confidence interval: −0.93 do −0.74: wiarygodny zakres prawdziwej korelacji. Często przydatniejszy niż p-value: nawet optymistyczny koniec tego przedziału to silna korelacja ujemna.
• sample estimates: cor = −0.87: ta sama liczba, którą dało cor().

Małe próby wymagają tu dodatkowej ostrożności: przy n = 10 korelacje rzędu ±0.5 pojawiają się przypadkiem niepokojąco często, a szeroki przedział ufności ci to pokaże. Podawaj przedział, nie tylko p-value.

Zobacz to: zawsze rysuj wykres

Współczynnik korelacji ściska całą zależność do jednej liczby, a to ściśnięcie może ukryć krzywiznę, skupiska albo jeden punkt, który robi całą robotę. Zanim zaufasz jakiemukolwiek r, obejrzyj wykres punktowy:

plot(mtcars$wt, mtcars$mpg)              # one pair
pairs(mtcars[, c("mpg", "wt", "hp", "disp")])  # every pair in the matrix

pairs() rysuje siatkę wykresów punktowych odpowiadającą macierzy korelacji, co jest najszybszym sposobem, żeby sprawdzić, czy liczby znaczą to, co myślisz. Do dopracowanych grafik macierzy w stylu mapy ciepła standardowym narzędziem jest pakiet corrplot (install.packages("corrplot"), potem corrplot(cor(m))).

Co warto zapamiętać

• cor(x, y) daje r Pearsona: znak to kierunek, wartość bezwzględna to siła, zawsze w [−1, 1], bez jednostek.
• Korelacja mierzy liniowe współzmienianie się i nie mówi nic o przyczynowości; ukryta trzecia zmienna zawsze jest kandydatem.
• method = "spearman" dla rang: dane porządkowe, wartości odstające, zależności monotoniczne, ale zakrzywione.
• cor(df) na kolumnach liczbowych daje macierz; zaokrąglij ją przez round(, 2) i pilnuj argumentu use =, gdy brakuje danych.
• cor.test(x, y) dodaje p-value i przedział ufności; podawaj przedział.
• Zawsze obejrzyj wykres punktowy, zanim uwierzysz liczbie.

Dalej: gdy pytanie zaostrza się z "czy zmieniają się razem?" do "czy średnia tej grupy różni się od tamtej?", czyli test t.