Корреляция в R: cor(), cor.test() и корреляционные матрицы

Корреляция между двумя переменными: cor()

Корреляция спрашивает: когда одна переменная растёт, вторая тоже склонна расти (положительная), падать (отрицательная) или жить своей жизнью (около нуля)? В R это один вызов:

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

Ответ — примерно −0.87: более тяжёлые машины расходуют больше топлива, и связь сильная. Это единственное число и есть коэффициент корреляции Пирсона, повсеместно обозначаемый r.

Как читать r

Коэффициент всегда лежит между −1 и +1. Знак даёт направление, величина даёт силу:

| |r| | Типичное прочтение | | --- | --- | | 0.0 – 0.2 | пренебрежимая | | 0.2 – 0.4 | слабая | | 0.4 – 0.6 | умеренная | | 0.6 – 0.8 | сильная | | 0.8 – 1.0 | очень сильная |

Относитесь к этим диапазонам как к поводу для разговора, а не как к закону: в физике r равное 0.6 разочаровывает, в психологии это событие карьеры. Два свойства стоит усвоить: у r нет единиц измерения (корреляция веса в тоннах с расходом даёт то же r, что и веса в килограммах, потому что r считается по стандартизованным значениям), и r измеряет только линейную связь — идеальная U-образная зависимость может иметь r ≈ 0.

И фраза, которую необходимо произнести: корреляция не есть причинность. Продажи мороженого и число утоплений сильно коррелируют по месяцам года — не потому, что мороженое топит людей, а потому, что лето двигает и то и другое. Корреляция говорит, что две переменные движутся вместе; о том, почему, она молчит. Может, x влияет на y, может, y влияет на x, а может, третья вещь (в данном случае сезон) влияет на обе. Решение между этими вариантами требует экспериментов или тщательного причинного рассуждения, а не большего r.

Спирмен и Кендалл: когда Пирсон — неподходящий инструмент

Пирсон работает по сырым значениям, что делает его чувствительным к выбросам и слепым к изогнутым связям. Аргумент method переключает на ранговые альтернативы:

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

Спирмен заменяет каждое значение его рангом, а затем считает Пирсона по рангам. Поскольку y всегда растёт с ростом x, все ранги выстраиваются, и Спирмен сообщает ровно 1 — размер выброса перестаёт иметь значение, важна только его позиция. Берите Спирмена, когда данные порядковые (шкалы опросов), сильно скошенные или когда связь монотонна, но не прямая. Кендалл отвечает на похожий вопрос через согласованные и несогласованные пары; он устойчивее на малых выборках, но медленнее, а вариантом по умолчанию обычно служит Спирмен.

Корреляционная матрица

Чтобы просканировать связи сразу по многим переменным, передайте cor() несколько числовых столбцов:

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

Каждая переменная против каждой, с единицами по диагонали (всё идеально коррелирует само с собой) и зеркальным отражением поперёк неё. Округление до двух знаков важнее, чем кажется: неокруглённая матрица это стена цифр, а смысл матрицы в том, чтобы её сканировать. Здесь взгляд показывает, что mpg отрицательно коррелирует со всеми тремя (более тяжёлые, мощные машины с большим двигателем сжигают больше топлива), а wt, hp и disp сильно положительно связаны между собой — кластер переменных «большой машины», который станет важен, когда вы дойдёте до линейной регрессии и её головной боли с мультиколлинеарностью.

Пропущенные значения: аргумент use

При пропущенных данных cor() по умолчанию возвращает NA, а не гадает:

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" отбрасывает каждую строку, содержащую любой NA, а затем считает всю матрицу по выжившим. Согласованно, но расточительно: пропущенный wt заодно удаляет эту строку и из пары mpg–hp.
• use = "pairwise.complete.obs" считает каждую ячейку по всем строкам, где присутствует именно эта пара. Так сохраняется больше данных, но разные ячейки опираются на разные подмножества, что изредка может дать матрицу, не согласованную внутри себя.

Для пары затесавшихся NA подойдёт любой вариант; просто укажите, каким вы воспользовались.

Значима ли она? cor.test()

cor() даёт число, но не даёт понимания, не является ли оно шумом. cor.test() проводит проверку гипотезы:

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

Пройдёмся по выводу блок за блоком:

• t = −9.56, df = 30 — статистика критерия. Нулевая гипотеза состоит в том, что истинная корреляция равна нулю; наблюдённое r переводится в статистику t с n − 2 степенями свободы (32 машины − 2).
• p-value = 1.29e-10 — если бы истинная корреляция была нулевой, вероятность увидеть r, настолько далёкое от нуля, в выборке из 32 наблюдений составляла бы около 0.0000000001. Это подавляющее свидетельство того, что связь реальна, — но помните, p-значение говорит о том, отличается ли r от нуля, а не о том, велика ли связь и причинна ли она.
• 95 percent confidence interval: −0.93 до −0.74 — правдоподобный диапазон для истинной корреляции. Часто он полезнее p-значения: даже оптимистичный конец этого интервала является сильной отрицательной корреляцией.
• sample estimates: cor = −0.87 — то же число, что дала cor().

Малые выборки заслуживают здесь особого уважения: при n = 10 корреляции ±0.5 возникают по случайности пугающе часто, и широкий доверительный интервал вам об этом скажет. Сообщайте интервал, а не только p-значение.

Увидеть своими глазами: всегда стройте график

Коэффициент корреляции сжимает целую связь в одно число, и это сжатие может скрыть изгиб, кластеры или одну точку, делающую всю работу. Прежде чем доверять любому r, посмотрите на диаграмму рассеяния:

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

pairs() рисует сетку диаграмм рассеяния, соответствующую вашей корреляционной матрице, — самый быстрый способ проверить, что числа означают то, что вы думаете. Для аккуратной графики матрицы в стиле тепловой карты стандартным инструментом служит пакет corrplot (install.packages("corrplot"), затем corrplot(cor(m))).

Что вы уносите с собой

• cor(x, y) даёт r Пирсона: знак это направление, величина это сила, всегда в диапазоне [−1, 1], без единиц измерения.
• Корреляция измеряет линейное совместное движение и ничего не говорит о причинности — притаившаяся третья переменная всегда является кандидатом.
• method = "spearman" для рангов: порядковые данные, выбросы, монотонные, но изогнутые связи.
• cor(df) по числовым столбцам даёт матрицу; округляйте её через round(, 2) и следите за аргументом use = при пропущенных данных.
• cor.test(x, y) добавляет p-значение и доверительный интервал — сообщайте интервал.
• Всегда смотрите на диаграмму рассеяния, прежде чем верить числу.

Дальше: когда вопрос заостряется от «движутся ли они вместе?» до «отличается ли среднее этой группы от среднего той?» — t-критерий.