14 Ковариация и корреляция

Возьмем новый набор данных, на этот раз про американских студентов, вес их рюкзаков и проблемы со спиной. Этот набор данных хранится в пакете Stat2Data — пакет с большим количеством разнообразных данных.

install.packages("Stat2Data")

С помощью функции data() загрузим набор данных Backpack:

library(tidyverse)
library(Stat2Data)
data(Backpack)

Давайте посмотрим, что внутри этой переменной:

skimr::skim(Backpack)
Таблица 14.1: Data summary
Name Backpack
Number of rows 100
Number of columns 9
_______________________
Column type frequency:
factor 3
numeric 6
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Major 0 1 FALSE 41 Bio: 9, Bus: 8, LS: 7, ME: 6
Sex 0 1 FALSE 2 Fem: 55, Mal: 45
Status 0 1 FALSE 2 U: 97, G: 3

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
BackpackWeight 0 1 11.66 5.77 2.00 8.00 11.00 14.25 35.00 ▅▇▂▁▁
BodyWeight 0 1 153.05 29.40 105.00 130.00 147.50 170.00 270.00 ▆▇▂▁▁
Ratio 0 1 0.08 0.04 0.02 0.05 0.07 0.10 0.18 ▆▇▆▃▁
BackProblems 0 1 0.32 0.47 0.00 0.00 0.00 1.00 1.00 ▇▁▁▁▃
Year 0 1 3.20 1.39 0.00 2.00 3.00 4.00 6.00 ▃▆▇▆▆
Units 0 1 14.27 2.81 0.00 13.00 15.00 16.00 19.00 ▁▁▁▇▆

С помощью ?Backpack можно получить подробное описание колонок этого датасета.

Например, можно заметить, что масса как самих студентов, так и их рюкзаков выражена в фунтах. Давайте создадим новые переменные backpack_kg и body_kg, в которых будет записан вес (рюкзаков и самих студентов соответственно) в понятным для нас килограммах. Новый набор данных сохраним под названием back.

back <- Backpack %>%
  mutate(backpack_kg = 0.45359237 * BackpackWeight,
         body_kg = 0.45359237 * BodyWeight)

До этого мы говорили о различиях между выборками. Теперь мы будем говорить о связи между переменными.

14.1 Ковариация

Самая простая мера связи между двумя переменными — это ковариация. Если ковариация положительная, то чем больше одна переменная, тем больше другая переменная. При отрицательной ковариации все наоборот: чем больше одна переменная, тем меньше другая.

  • Формула ковариации:

\[\sigma_{xy} = cov(x, y) = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{n}\]

  • Оценка ковариации по выборке:

\[\hat{\sigma}_{xy} = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{n-1}\]

Ковариация переменной самой с собой — дисперсия.

В R есть функция cov() для подсчета ковариации. На самом деле, функция var() делает то же самое. Обе эти функции считают сразу матрицу ковариаций для всех сочетаний колонок на входе:

back %>%
  select(body_kg, backpack_kg) %>%
  cov()
##                body_kg backpack_kg
## body_kg     177.807700    6.601954
## backpack_kg   6.601954    6.838333
back %>%
  select(body_kg, backpack_kg) %>%
  var()
##                body_kg backpack_kg
## body_kg     177.807700    6.601954
## backpack_kg   6.601954    6.838333

Ну а по углам этой матрицы — дисперсии!

14.2 Корреляция

Корреляцией обычно называют любую связь между двумя переменными, это просто синоним слова “ассоциация.” Если вдруг слово “корреляция” вам еще непривычно, то попробуйте мысленно заменять “корреляцию” на “ассоциацию,” а “коррелирует” на “связано.” Коэффициент корреляции — это уже конкретная математическая формула, которая позволяет посчитать эту связь и принимает значения от -1 до 1.37

  • Если коэффициент корреляции положительный, то чем больше значения в одной переменной, тем больше значения в другой переменной.

  • Если коэффициент корреляции отрицательный, то чем больше значения в одной переменной, тем меньше значения в другой переменной.

  • Если коэффициент корреляции равен 0, то изменения одной переменной не связано с изменениями в другой переменной.

14.2.1 Коэффициент корреляции Пирсона

Самый известный коэффициент корреляции - коэффициент корреляции Пирсона:

\[\rho_{xy} = \frac{\sigma_{xy}}{\sigma_x \sigma_y} = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{\sqrt{\sum_{i = 1}^n(x_i - \overline{x})^2}\sqrt{\sum_{i = 1}^n(y_i - \overline{y})^2}} = \frac{1}{n}\sum_{i = 1}^n z_{x,i} z_{y, i}\]

Оценка коэффициента корреляции Пирсона по выборке: \[r_{xy} = \frac{\hat{\sigma}_{xy}}{\hat{\sigma}_x \hat{\sigma}_y} = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{\sqrt{\sum_{i = 1}^n(x_i - \overline{x})^2}\sqrt{\sum_{i = 1}^n(y_i - \overline{y})^2}} = \frac{1}{n - 1}\sum_{i = 1}^n z_{x,i} z_{y, i}\]

Коэффициент корреляции Пирсона можно понимать по-разному. С одной стороны, это просто ковариация, нормированная на стандартное отклонение обоих переменных. С другой стороны, можно понимать это как среднее произведение z-оценок.

Корреляцию в R можно посчитать с помощью функции cor():

back %>%
  select(body_kg, backpack_kg) %>%
  cor()
##               body_kg backpack_kg
## body_kg     1.0000000   0.1893312
## backpack_kg 0.1893312   1.0000000

Для тестирования уровня значимости нулевой гипотезы для корреляции есть функция cor.test(). В случае с коэффициентами корреляции, нулевая гипотеза формулируется как отсутствие корреляции (т.е. она равна нулю) в генеральной совокупности.

cor.test(back$backpack_kg, back$body_kg)
## 
##  Pearson's product-moment correlation
## 
## data:  back$backpack_kg and back$body_kg
## t = 1.9088, df = 98, p-value = 0.05921
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.007360697  0.371918344
## sample estimates:
##       cor 
## 0.1893312

Результат выполнения этой функции очень похож на то, что мы получали при проведении t-теста.

14.2.2 Непараметрические коэффициенты корреляции

У коэффициента корреляции Пирсона, как и у t-теста, есть свои непараметрические братья: коэффициент корреляции Спирмена и коэффициент корреляции Кэнделла. Из них чаще используется коэффициент корреляции Спирмена. Посчитать его можно с помощью той же функции cor.test(), задав соответствующее значение параметра method =:

cor.test(back$backpack_kg, back$body_kg, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  back$backpack_kg and back$body_kg
## S = 131520, p-value = 0.03527
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2108001
cor.test(back$backpack_kg, back$body_kg, method = "kendall")
## 
##  Kendall's rank correlation tau
## 
## data:  back$backpack_kg and back$body_kg
## z = 2.083, p-value = 0.03725
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.1478736

Заметьте, в данном случае два метода хотя и привели к схожим размерам корреляции, но в одном случае p-value оказался больше 0.05, а в другом случае - меньше 0.05. Выбирать тест a posteriori на основе того, какие результаты вам нравятся больше, — плохая практика (@ref(bad_practice)). Не надо так делать.

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

Возможно, вы нашли что-то более интересное для проверки гипотезы о корреляции. Например, вы еще хотите проверить гипотезу о связи количества учебных кредитов и массе рюкзака: логично предположить, что чем больше студент набрал себе курсов, тем тяжелее его рюкзак (из-за большего количества учебников). Или что студенты к старшим курсам худеют и становятся меньше. Или что те, кто набрал себе много курсов, меньше питаются и от того меньше весят. В общем, хотелось бы прокоррелировать все интересующие нас переменные со всеми. Это можно сделать с помощью функции cor():

back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  cor()
##                 body_kg backpack_kg       Units        Year
## body_kg      1.00000000  0.18933115 -0.23524088 -0.09301727
## backpack_kg  0.18933115  1.00000000  0.09438453  0.05762194
## Units       -0.23524088  0.09438453  1.00000000 -0.02946373
## Year        -0.09301727  0.05762194 -0.02946373  1.00000000

Но функция cor()не позволяет посчитать p-value для этих корреляций! Функция cor.test() позволяет получить p-value, но только для одной пары переменных.

На помощь приходит пакет psych с функцией corr.test():

back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  psych::corr.test()
## Call:psych::corr.test(x = .)
## Correlation matrix 
##             body_kg backpack_kg Units  Year
## body_kg        1.00        0.19 -0.24 -0.09
## backpack_kg    0.19        1.00  0.09  0.06
## Units         -0.24        0.09  1.00 -0.03
## Year          -0.09        0.06 -0.03  1.00
## Sample Size 
## [1] 100
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##             body_kg backpack_kg Units Year
## body_kg        0.00        0.30  0.11    1
## backpack_kg    0.06        0.00  1.00    1
## Units          0.02        0.35  0.00    1
## Year           0.36        0.57  0.77    0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Тем не менее, если у вас много гипотез для тестирования, то у вас появляется проблема: вероятность выпадения статистически значимых результатов сильно повышается. Даже если эти переменные никак не связаны друг с другом.

Эта проблема называется проблемой множественных сравнений (multiple comparisons problem).38 Если мы проверяем сразу несколько гипотез, то у нас возрастает групповая вероятность ошибки первого рода (Family-wise error rate) — вероятность ошибки первого рода для хотя бы одной из множества гипотез.

Например, если вы коррелируете 10 переменных друг с другом, то вы проверяете 45 гипотез о связи. Пять процентов из этих гипотез, т.е. в среднем 2-3 гипотезы у вас будут статистически значимыми даже если никаких эффектов на самом деле нет!

Поэтому если вы проверяете сразу много гипотез, то необходимо применять поправки на множественные сравнения (multiple testing correction). Эти поправки позволяют контролировать групповую вероятность ошибки первого рода на желаемом уровне. Самая простая и популярная поправка на множественные сравнения — поправка Бонферрони (Bonferroni correction). Она считается очень просто: мы просто умножаем p-value на количество проверяемых гипотез!

back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  psych::corr.test(adjust = "bonferroni")
## Call:psych::corr.test(x = ., adjust = "bonferroni")
## Correlation matrix 
##             body_kg backpack_kg Units  Year
## body_kg        1.00        0.19 -0.24 -0.09
## backpack_kg    0.19        1.00  0.09  0.06
## Units         -0.24        0.09  1.00 -0.03
## Year          -0.09        0.06 -0.03  1.00
## Sample Size 
## [1] 100
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##             body_kg backpack_kg Units Year
## body_kg        0.00        0.36  0.11    1
## backpack_kg    0.06        0.00  1.00    1
## Units          0.02        0.35  0.00    1
## Year           0.36        0.57  0.77    0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Это очень “дубовая” и излишне консервативная поправка. Да, она гарантирует контроль групповую вероятности ошибки первого рода, но при этом сильно повышает вероятность ошибки второго рода — вероятность пропустить эффект, если он на самом деле существует. Поэтому по умолчанию в R используется более либеральная поправка на множественные сравнения под названием поправка Холма или поправка Холма-Бонферрони (Holm-Bonferroni correction), которая, тем не менее, тоже гарантирует контроль групповой вероятности ошибки первого рода.

Альтернативный подход к решению проблемы множественных сравнений — это контроль средней доли ложных отклонений (False Discovery Rate; FDR) на на уровне не выше уровня \(\alpha\). Это более либеральный подход: в данном случае мы контролируем, что ложно-положительных результатов у нас не больше, например, 5%. Такой подход применяется в областях, где происходит масштабное множественное тестирование. Попытка контролировать групповую вероятность ошибки первого уровня не выше уровня \(\alpha\) привела бы к чрезвычайно низкой вероятности обнаружить хоть какие-нибудь эффекты (т.е. к низкой статистической мощности).

Самая известная поправка для контроля средней доли ложных отклонений — это поправка Бенджамини — Хохберга (Benjamini-Hochberg correction).

back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  psych::corr.test(adjust = "BH")
## Call:psych::corr.test(x = ., adjust = "BH")
## Correlation matrix 
##             body_kg backpack_kg Units  Year
## body_kg        1.00        0.19 -0.24 -0.09
## backpack_kg    0.19        1.00  0.09  0.06
## Units         -0.24        0.09  1.00 -0.03
## Year          -0.09        0.06 -0.03  1.00
## Sample Size 
## [1] 100
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##             body_kg backpack_kg Units Year
## body_kg        0.00        0.18  0.11 0.54
## backpack_kg    0.06        0.00  0.54 0.68
## Units          0.02        0.35  0.00 0.77
## Year           0.36        0.57  0.77 0.00
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Все перечсиленные поправки (и еще несколько других) доступны не только в функции corr.test(), но и в базовом R с помощью функции p.adjust(). Эта функция принимает вектор из p-value и возвращает результат применения поправок.

p_vec <- seq(0.0001, 0.06, length.out = 10)
p_vec
##  [1] 0.000100000 0.006755556 0.013411111 0.020066667 0.026722222 0.033377778
##  [7] 0.040033333 0.046688889 0.053344444 0.060000000
p.adjust(p_vec) #по умолчанию используется поправка Холма-Бонферрони
##  [1] 0.0010000 0.0608000 0.1072889 0.1404667 0.1603333 0.1668889 0.1668889
##  [8] 0.1668889 0.1668889 0.1668889
p.adjust(p_vec, method = "bonferroni")
##  [1] 0.00100000 0.06755556 0.13411111 0.20066667 0.26722222 0.33377778
##  [7] 0.40033333 0.46688889 0.53344444 0.60000000
p.adjust(p_vec, method = "BH")
##  [1] 0.00100000 0.03377778 0.04470370 0.05016667 0.05344444 0.05562963
##  [7] 0.05719048 0.05836111 0.05927160 0.06000000

14.4 Хитмэп корреляций

  1. Впрочем, это не так уж и важно: на практике часто опускают слово “коэффициент” и называют корреляцией как и просто связь, так и ее способ измерения.↩︎

  2. “Проблема множественных сравнений” - это устоявшийся термин, который используется и в случае множественных корреляций, и в случае множественных сравнений средних и в любых других случаях с тестированием нескольких гипотез одновременно.↩︎