330-cov_cor.qmd

# Ковариация и корреляция {#sec-cov_cor}

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

```{r, eval = FALSE}
install.packages("Stat2Data")
```

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

```{r}
library(tidyverse)
library(Stat2Data)
data(Backpack)
```

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

```{r}
skimr::skim(Backpack)
```

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

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

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

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

## Ковариация {#sec-cov}

Самая простая мера связи между двумя переменными --- это **ковариация *(covariation).*** Если **ковариация положительная**, то чем *больше* одна переменная, тем *больше* другая переменная. При **отрицательной ковариации** все наоборот: чем *больше* одна переменная, тем *меньше* другая. **Ковариация** между переменными $x$ и $y$ часто обозначается буквой $\sigma_{xy}$. Похоже на дисперсию и стандартное отклонение (см. @sec-var), не правда ли? И скоро мы убедимся, что это не случайно!

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

$$\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()` для подсчета ковариации. Эта функция считает сразу **матрицу ковариаций *(covariation matrix)*** для всех сочетаний колонок на входе:

```{r}
back %>%
  select(body_kg, backpack_kg) %>%
  cov()
```

Полученная матрица симметрична: порядок двух переменных не имеет значения. Но что мы имеем по основной диагонали -- это, получается, ковариация переменной с самой собой? Давайте подставим в формулу для ковариации $x$ вместо $y$, чтобы узнать, что скрывается под ковариацией переменной с самой собой:

$$\sigma_{xx} = cov(x, x) = \frac{\sum_{i = 1}^n(x_i - \overline{x})(x_i - \overline{x})}{n} = \frac{\sum_{i = 1}^n(x_i - \overline{x})^2}{n}$$

Знакомая формула? Это формула дисперсии (см. @sec-var)! Выходит, что **ковариация переменной с самой собой -- это дисперсия.** Таким образом, **в матрице ковариаций по основным диагоналям находятся дисперсии переменных.**

$$
\begin{bmatrix}\sigma_x^2 & \sigma_{xy}\\\sigma_{xy} & \sigma_y^2\end{bmatrix}
$$

::: callout-tip
## *Полезное:* variance и covariance

В английском языке **дисперсия --** ***variance***, а **ковариация --** ***covariance***, что значительно упрощает запоминание связи этих двух понятий.
:::

Что интересно, в R это тоже отображено: посчитать матрицу ковариаций можно как с помощью функции `cov()`, так и уже знакомой нам функции `var()` для расчета дисперсии. Различаются эти функции только дополнительными параметрами.

```{r}
back %>%
  select(body_kg, backpack_kg) %>%
  cov()

back %>%
  select(body_kg, backpack_kg) %>%
  var()
```

::: callout-warning
## *Для продвинутых:* Ковариационная матрица с точки зрения линейной алгебры

Расчет матрицы ковариаций легко представить в виде матричных операций:

1.  Сначала мы вычтем из каждого столбца среднее по столбцу:

```{r}
X <- back %>%
  select(body_kg, backpack_kg) %>%
  as.matrix() %>%
  apply(2, function(x) x - mean(x))
```

2.  Затем транспонируем матрицу (поменяем местами столбцы и строки) и перемножим саму на себя: $$
    X^TX
    $$

```{r}
t(X) %*% X
```

3.  Поделим матрицу на $n$ или $n - 1$ (если мы хотим оценить матрицу ковариаций в генеральной совокупности по выборке), где $n$ -- количество наблюдений. Это и будет ковариационной матрицей!

$$
C_X = \frac{X^TX}{n-1}
$$

```{r}
t(X) %*% X / (nrow(back) - 1)
back %>%
  select(body_kg, backpack_kg) %>%
  cov()
```
:::

Ковариации имеют большое значение в статистике. Расчет ковариационной матрицы -- это необходимый этап для многих многомерных методов, например, для анализа главных компонент (см. @sec-dim_red_pca). Однако у ковариации есть серьезное ограничение -- ее размер привязан к исходной шкале, поэтому сложно оценить, насколько ковариация большая или маленькая. Поэтому на практике гораздо больше используются коэффициенты корреляции.

## Корреляция {#sec-cor}

**Корреляцией** обычно называют любую связь между двумя переменными, это просто синоним слова "ассоциация". Если вдруг слово "корреляция" вам еще непривычно, то попробуйте мысленно заменять "корреляцию" на "ассоциацию", а "коррелирует" на "связано". **Коэффициент корреляции** --- это уже конкретная математическая формула, которая позволяет посчитать эту связь и принимает значения от $-1$ до $1$.[^330-cov_cor-1]

[^330-cov_cor-1]: Впрочем, это не так уж и важно: на практике часто опускают слово "коэффициент" и называют корреляцией как и просто связь, так и ее способ измерения.

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

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

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

### Коэффициент корреляции Пирсона {#sec-pearson_cor}

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

$$\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$-оценок (@sec-z_scores).

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

```{r}
back %>%
  select(body_kg, backpack_kg) %>%
  cor()
```

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

```{r}
cor.test(back$backpack_kg, back$body_kg)
```

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

### Непараметрические коэффициенты корреляции {#sec-nonparam_cor}

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

```{r}
cor.test(back$backpack_kg, back$body_kg, method = "spearman")
cor.test(back$backpack_kg, back$body_kg, method = "kendall")
```

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

## Корреляционная матрица {#sec-cor_mat}

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

```{r}
back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  cor()
```

В корреляционной матрице все значения находятся в диапазоне от $-1$ до $1$, а по основной диагонали находятся единицы:

$$
\begin{bmatrix}1 & \rho_{xy}\\\rho_{xy} & 1\end{bmatrix}
$$

И действительно: переменная коррелирует сама с собой идеально, то есть коэффициент корреляции равен $1$.

::: callout-warning
## *Для продвинутых:* Корреляционная матрица с точки зрения линейной алгебры

Корреляционную матрицу можно построить разными способами, но самый простой -- сделать $z$-преобразование переменных (@sec-z_scores), а для полученных z-оценок посчитать ковариационную матрицу.

1.  Произведем $z$-преобразования каждого столбца:

```{r}
X <- back %>%
  select(body_kg, backpack_kg) %>%
  as.matrix() %>%
  scale()
```

2.  Затем транспонируем полученную матрицу и перемножим саму на себя по правилам линейной алгебры: $$
    X^TX
    $$

```{r}
t(X) %*% X
```

3.  Осталось поделить результаты матрицу на $n - 1$, где $n$ -- количество наблюдений. Это и будет ковариационной матрицей!

$$
COR_X = \frac{X^TX}{n-1}
$$

```{r}
t(X) %*% X / (nrow(back) - 1)
back %>%
  select(body_kg, backpack_kg) %>%
  cov()
```
:::

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

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

```{r}
back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  psych::corr.test()
```

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

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

[^330-cov_cor-2]: "Проблема множественных сравнений" - это устоявшийся термин, который используется и в случае множественных корреляций, и в случае множественных сравнений средних и в любых других случаях с тестированием нескольких гипотез одновременно.

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

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

```{r}
back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  psych::corr.test(adjust = "bonferroni")
```

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

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

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

```{r}
back %>%
  select(body_kg, backpack_kg, Units, Year) %>%
  psych::corr.test(adjust = "BH")
```

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

```{r}
p_vec <- seq(0.0001, 0.06, length.out = 10)
p_vec
p.adjust(p_vec) #по умолчанию используется поправка Холма-Бонферрони
p.adjust(p_vec, method = "bonferroni")
p.adjust(p_vec, method = "BH")
```

## Хитмэп корреляций {#sec-hitmap_cor}

Как видите, почти все коррелирует друг с другом, даже с учетом поправок. Такие множественные корреляции лучше всего смотреть с помощью хитмап-визуализации.

В качестве примера возьмем встроенный датасет `mtcars`, в котором есть множество количественных переменных.

```{r}
mtcars
```

Для визуализации возьмем пакет `{corrplot}`.

```{r, eval = FALSE}
install.packages("corrplot")
```

Для начала нам нужно построить матрицу корреляций. Уже знакомая нам функция `cor()` для это вполне подойдет:

```{r}
library(corrplot)
cor(mtcars)
```

Основная функция пакета `{corrplot}` -- функция `corrplot()`. Давайте теперь попробуем ее применить на нашей матрице корреляций.

```{r}
corrplot(cor(mtcars))
```

По умолчанию значение корреляций кодируется цветом и размером круга. У функции `corrplot()` есть множество параметров, позволяющих довольно тонко настраивать хитмэп корреляций. Например, можно кодировать только цветом, а переменные сгруппировать на основе кластеризации.

```{r}
corrplot(cor(mtcars), method = "color", order = "hclust")
```

Можно сделать еще интереснее: добавить крестики, которые будут означать наличие или отсутствие статистической значимости. В этом случае нам понадобится не только матрица корреляций, но соответствующая матрица с *p-value*. Для этого нам нужно обратиться к уже знакомой нам функции `corr.test()` из пакета `{psych}`.

```{r}
mtcars_cor_p <- psych::corr.test(mtcars)
mtcars_cor_p
```

Нам нужно разворошить полученную переменную mtcars_cor_p

```{r}
str(mtcars_cor_p)
```

```{r}
corrplot(corr = mtcars_cor_p$r,
         p.mat = mtcars_cor_p$p,
         method = "color",
         order = "hclust")
```

Крестиками отмечены статистические незначимые (с *p.value \< .05),* а без крестиков -- статистически значимые коэффициенты корреляции. Обратите внимание, матрица несимметричная: `psych::corr.test()` возвращает скорректированные *p-value* в верхнем правом треугольнике в матрице корреляций.

## Матрица корреляций в виде графа

```{r, eval = FALSE}
install.packages("corrr")
```

```{r}
library(corrr)
```

Функционал пакета {corrr} позволяет работать с корреляционными матрицами в духе tidyverse, возвращая тиббл вместо матрицы.

```{r}
correlate(mtcars)
```

Нас же интересует в первую очередь представление матрицы корреляций в виде графа.

Делается это следующим образом: мы устанавливаем минимальное значение для корреляций с помощью параметра `min_cor =`, которые будут отображены в качестве связи между переменными. Если мы этого не сделаем, то мы просто получим полностью связанный граф, потому что все переменные хотя бы немного коррелируют со всеми. По умолчанию этот параметр равен 0.3, но мы можем изменить это значение на собственное усмотрение. В данном случае установим `min_cor = 0.7`, потому что в нашей матрице очень много сильных корреляций.

```{r}
correlate(mtcars) %>%
  corrr::network_plot(min_cor = 0.7)
```