Коэффициенты корреляции

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

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

Применение

  • используется, когда оба показателя измерены в количественной шкале (примеры: рост в сантимерах и вес в килограммах, число лет опыта работы и заработная плата);
  • чувствителен к выбросам, нехарактерным значениям, которые можно увидеть на диаграмме рассеяния.

Расчет коэффициента

\[ r = \frac{\sum_{i=1}^n(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2 \sum_{i=1}^n(y_i-\bar{y})^2}}, \] где \(n\) — число наблюдений в выборках (число пар наблюдений), \(\bar{x}\) — среднее арифметическое по выборке \(X\), \(\bar{y}\) — среднее арифметическое по выборке \(Y\).

Маленький пример (4 наблюдения — это очень мало, просто чтобы быстро посчитать):

##   x  y
## 1 1  2
## 2 4  2
## 3 5  9
## 4 8 19

Посчитаем средние (\(\bar{x}\) и \(\bar{y}\)):

mean(x); mean(y)
## [1] 4.5
## [1] 8

Вычтем из каждого значения \(x\) среднее:

x - mean(x)
## [1] -3.5 -0.5  0.5  3.5

Проделаем то же самое для \(y\):

y - mean(y)
## [1] -6 -6  1 11

Перемножим результаты:

(x - mean(x)) * (y - mean(y))
## [1] 21.0  3.0  0.5 38.5

Сложим — числитель дроби готов!

sum((x - mean(x)) * (y - mean(y)))
## [1] 63

Посчитаем знаменатель:

sum((x - mean(x))^2) # первая сумма
## [1] 25
sum((y - mean(y))^2) # вторая сумма
## [1] 194
sqrt(sum((x - mean(x))^2) * sum((y - mean(y))^2)) # итог
## [1] 69.64194

Посчитаем коэффициент:

63 / 69.64194   # r
## [1] 0.9046273

Коэффициент положительный, связь между показателями прямая (чем больше \(x\), тем больше \(y\)). Коэффициент высокий, по модулю близок к 1, связь сильная.

Нулевая и альтернативная гипотезы

\[ H_0: \rho = 0 \text{ (связи нет)} \] \[ H_1: \rho \ne 0 \text{ (связь есть)} \]

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

Значение \(\rho\) в гипотезах — значение коэффициента корреляции между двумя генеральными совокупностями, не конкретными выборками.

Если отвергаем гипотезу, то можем заключить, что связь между показателями есть.

Расчет коэффициента и проверка его значимости в R

cor.test(d$x, d$y)
## 
##  Pearson's product-moment correlation
## 
## data:  d$x and d$y
## t = 3.0017, df = 2, p-value = 0.09537
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4323956  0.9980148
## sample estimates:
##       cor 
## 0.9046273

Согласно выдаче, коэффициент корреляции Пирсона равен 0.905 (значение cor в sample estimates). Значение p-value достаточно велико, больше 0.05, поэтому у нас нет оснований отвергнуть нулевую гипотезу. Истинное значение коэффициента корреляции равно 0, связи между показателями нет, коэффициент незначим. В нашем случае это объяснимо: в датасете всего 4 наблюдения, на таком объеме сложно получить значимые результаты.

Что еще есть в выдаче?

  • t: значение статистики критерия, которая используется для проверки гипотезы о равенстве коэффициента корреляции нулю. Имеет t-распределение с числом степеней свободы \(df = n - 2 = 4 - 2 = 2\) (в выдаче это df). Считается так:

\[ t = r \times \sqrt{\frac{n-2}{1 - r^2}} = 0.905 \times \sqrt{\frac{4-2}{1-0.905^2}} \approx 3.008. \]

  • 95 percent confidence interval: 95% доверительный интервал для коэффициента корреляции.

Коэффициент корреляции Ч.Спирмена

Применение

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

Расчет коэффициента

\[ r_{\text{Спирмена}} = 1 - \frac{6 \times \sum_{i=1}^n d_i^2}{n(n^2 -1)}, \]

где \(d_i^2\) — сумма квадратов разностей рангов наблюдений, \(n\) — число наблюдений в выборках (число пар наблюдений).

Маленький пример (4 наблюдения — это очень мало, просто чтобы быстро посчитать):

##   x  y
## 1 1  2
## 2 4  2
## 3 5  9
## 4 8 19

Посчитаем ранги наблюдений в x и y:

rank(x)
## [1] 1 2 3 4
rank(y)
## [1] 1.5 1.5 3.0 4.0

Посчитаем их разности и возведем их в квадрат:

(rank(x) - rank(y))^2 
## [1] 0.25 0.25 0.00 0.00

Суммируем квадраты разностей рангов:

sum((rank(x) - rank(y))^2) 
## [1] 0.5

Число пар наблюдений равно \(n = 4\) (четыре строки в датасете), подставим все в формулу и рассчитаем коэффициент:

1 - 6 * 0.5 / (4 * (4^2 - 1))
## [1] 0.95

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

Нулевая и альтернативная гипотезы

\[ H_0: \text{признаки независимы (связи нет)} \]

\[ H_1: \text{признаки не независимы (связи нет)} \] Приведенная выше нулевая гипотеза также называется гипотезой о значимости коэффициента корреляции.

Расчет коэффициента и проверка его значимости в R

cor.test(d$x, d$y, method = 'spearman')  # можно сократить до 'sp'
## Warning in cor.test.default(d$x, d$y, method = "spearman"): Cannot compute
## exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  d$x and d$y
## S = 0.51317, p-value = 0.05132
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9486833

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

Согласно выдаче, коэффициент корреляции Спирмена равен 0.948 (значение cor в sample estimates). Значение p-value чуть больше 0.05, поэтому у нас нет оснований отвергнуть нулевую гипотезу (хотя тут ситуация спорная, зависисит от исследователя). Признаки независимы, коэффициент незначим.

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

Парная (простая линейная регрессия)

Постановка задачи

Возьмем данные из файла english.csv, который использовался в практикуме и построим модель линейной взаимосвязи (линейную регрессию), которая бы объясняла, каким образом время, затраченное на угадывание слова (RTlexdec) зависит от его встречаемости в письменных текстах (WrittenFrequency).

english <- read.csv("http://math-info.hse.ru/f/2018-19/ling-data/english.csv")

Итак:

  • зависимая переменная: RTlexdec;
  • независимая переменная: WrittenFrequency.

Примечание 1. Обычно, исходя из задачи ясно, какая переменная является зависимой, а какая — независимой. Например, в нашем случае ясно, что встречаемость слова влияет на скорость его узнавания, не наоборот: вряд ли тот факт, что участники некоторого эксперимента быстро узнают слово, сделает это слово более популярным в текстах на английском языке. Еще пример: в паре переменных «доход человека» и «среднемесячные траты на обеды в кафе» первая является независимой, а вторая — зависимой, потому что очевидно, что более богатый человек может позволить себе тратить больше денег на обеды вне дома, но от того, что человек станет тратить больше денег на еду, богаче он не станет.

Примечание 2: Хотя мы и говорим «зависимая и независимая переменные», «влияние одной переменной на другую», регрессия, как и корреляция не показывает причинно-следственную связь, а просто позволяет выяснить характер связи между показателями.

Линейная модель

Уравнение модели (уравнение прямой):

\[ \text{RTlexdec} = \beta_0 + \beta_1 \times \text{WrittenFrequency} \] В отличии от корреляции, здесь мы рассматриваем направленную связь, то есть у нас есть предположение о том, какая из двух переменных оказывает влияние на другую. При этом мы хотим не просто выяснить, каким образом две переменные связаны (прямой или обратной связью), но еще хотим понять, насколько, в целом, изменяется значение зависимой переменной, когда независимая изменяется на единицу. Ответ на этот вопрос нам даст значение \(\beta_1\), которое будет рассчитано в R. Интерпретация уравнения выше:

  • При увеличении WrittenFrequency на единицу, значение RTlexdec, в среднем, изменяется на \(\beta_1\) (уменьшается или увеличивается — зависит от знака коэффициента.

  • Среднее (среднее ожидаемое) значение RTlexdec равно \(\beta_1\).

Проверка гипотез

Помимо значения коэффициента \(\beta_1\) (константа \(\beta_0\) нас обычно не интересует), необходимо знать, является ли этот коэффициент статистически значимым, то есть правда, что истинный коэффициент регрессии для генеральной совокупности будет существенно отличен от нуля. Если коэффициент статистически значим, то связь между переменными действительно есть, коэффицент имеет смысл принимать в рассмотрение и интерпретировать.

Нулевая гипотеза (она же гипотеза о значимости коэффициента регрессии):

\[ H_0: \beta_1 = 0 \text{ (независимая переменная не влияет на зависимую)} \]

Альтернативная гипотеза:

\[ H_1: \beta_1 \ne 0 \text{ (независимая переменная влияет на зависимую)} \]

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

Построение модели в R

Построим модель, описанную выше, используя функцию lm() — от linear model. На первом месте в уравнении всегда идет зависимая переменная, на втором, после ~ (тильда) — независимая. Далее запросим summary() по модели, чтобы узнать про нее все.

model <- lm(data = english, RTlexdec ~ WrittenFrequency)
summary(model)
## 
## Call:
## lm(formula = RTlexdec ~ WrittenFrequency, data = english)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -305.21  -83.71   -8.39   69.81  559.05 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      843.5759     4.4456  189.76   <2e-16 ***
## WrittenFrequency -26.9744     0.8311  -32.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 103.5 on 4566 degrees of freedom
## Multiple R-squared:  0.1874, Adjusted R-squared:  0.1873 
## F-statistic:  1053 on 1 and 4566 DF,  p-value: < 2.2e-16

Теперь обо всем по порядку.

1. Строка Call напоминает нам формулу для построенной модели.

2. Строка Residuals содержит описательные статистики для остатков (ошибок) модели: \[ \varepsilon = y - \hat{y}, \] где \(y\) — наблюдаемые значения зависимой переменной (из датасета), а \(\hat{y}\) — значения зависимой переменной, предсказанные моделью.

3. Таблица Coefficients — самая главная. В ней содержатся коэффициенты модели (\(\beta_0\) и \(\beta_1\)), а также p-value по итогам проверки значимости коэффициентов.

Intercept в первом столбце означает, что первая строка в таблице содержит характеристики константы (\(\beta_0\)), а WrittenFrequency во второй — что здесь содержатся характеристики коэффициента при этой переменной (\(\beta_1\)). В столбце Estimate находятся значения этих коэффициентов. В столбце Pr(>|t|) — соответствующие p-value из проверки гипотезы о значимости коэффициента регрессии.

Запишем уравнение данной модели на основе выдачи:

\[ \text{RTlexdec} = 843.58 - 26.97 \times \text{WrittenFrequency}. \] Если мы посмотрим на p-value для коэффициента при WrittenFrequency, мы увидим, что оно очень мало, следовательно, гипотезу о равенстве истинного коэффициента регресии нулю нужно отвергнуть. Коэффициент при WrittenFrequency значим, отличен от нуля. Поэтому его можно проинтерпретировать:

  • При увеличении встречаемости слова в текстах (WrittenFrequency) на единицу, время узнавания слова (RTlexdec) в среднем уменьшается на 26.97 миллисекунд.

4. Еще в таблице Coefficients есть два столбца: Std. Error и t value. Столбец Std. Error содержит значения стандартных ошибок коэффициентов, то есть оценок стандартных отклонений коэффициентов, полученных на основе наших данных. В столбце записаны значения статистик критерия, используемого для проверки гипотезы о равенстве коэффициента регрессии нулю. Статистики имеют t-распределение, распределение Стьюдента с числом степеней свободы \(\text{df} = n - 1\). Считаются они просто:

\[ t_0 = \frac{\beta_0}{se(\beta_0)} \text{ и } t_1 = \frac{\beta_1}{se(\beta_1)}, \]

где \(se(\beta_0)\) и \(\beta_1\) — стандартные ошибки коэффициентов.

Проверим на нашей выдаче:

\[ t_0 = \frac{843.5759}{4.4456} = 189.76 \text{ и } t_0 = \frac{-26.9744}{0.8311} = -32.45. \]

Строка с Signif. codes представляет собой расшифровку обозначений. Звездочки отвечают за значимость коэффициентов (как и p-value). Если

  • p-value менее 0.001 (***), то коэффициент значим на уровне значимости 0.1%;
  • p-value менее 0.01 (**), то коэффициент значим на уровне значимости 1%;
  • p-value менее 0.05 (*), то коэффициент значим на уровне значимости 5%;
  • p-value менее 0.1 (.), то коэффициент значим на уровне значимости 10%.

Проверка качества модели

1. Одним из показателей качества модели считается коэффициент детерминации \(R^2\). Он обозначает долю дисперсии зависимой переменной, которая объясняется моделью. Выражаясь неформально, \(R^2\) показывает, какую долю «реальности» объясняет модель. Чем ближе этот показатель к 1, тем лучше. Однако не существует единого соглашения о том, какой \(R^2\) считать высоким, зависит от исследовательской области и от данных. Для определенности можем считать, что значения \(R^2\) выше 0.6 являются высокими.

В парной регрессии \(R^2\) считается как коэффициент Пирсона между зависимой и независимой переменной, возведенный в квадрат:

\[ R^2 = r^2. \] В выдаче R коэффициент детерминации обозначен Multiple R-squared. В нашей модели \(R^2\) невысокий.

2. Чтобы оценки коэффициентов, полученные с помощью метода наименьших квадратов (а именно он и реализуется по умолчанию) были точными, должен выполняться ряд условий или допущений модели. Эти условия называют условиями Гаусса-Маркова.

Условия Гаусса-Маркова

  1. линейная зависимость между зависимой и независимой переменной;
  2. отсутствие закономерностей (паттернов) в распределении остатков модели, они должны быть случайны;
  3. математическое ожидание остатков равно нулю, \(E(\varepsilon) = 0\);
  4. дисперсия остатков постоянна, не изменяется при изменении значений независимой переменной, \(D(\varepsilon) = \sigma^2\);
  5. распределение остатков модели является нормальным.

Последнее условие обычно не относят к условиям Гаусса-Маркова. Если оно не выполняется, то оценки коэффиентов, полученные методом наименьших квадратов, все равно считаются лучшими среди несмещенных линейных оценок (англ. BLUE — Best linear unbiased estimators), но если это условие выполняется, то оценки коэффициентов становятся лучшими не только среди линейных, но среди всех (англ. BUE — Best Unbiased Estimators). Но так как в нашем курсе мы рассматриваем гауссовскую модель, где нормальность остатков предполагается, давайте будем учитывать последнее условие.

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

library(ggplot2)
ggplot(data = english, aes(x = WrittenFrequency, y = RTlexdec)) +
  geom_point()

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

ggplot(data = english, aes(x = WrittenFrequency, y = RTlexdec)) +
  geom_point() + geom_smooth(method = lm)

Второе, третье и четвертое условия тоже легко проверить графически — нужно построить диаграмму рассеяния «остатки vs независимая переменная». А для этого нужно вытащить значения остатков из модели и добавить их отдельным столбцом в датасет:

english$residuals <- model$residuals

Теперь строим диаграмму рассеяния и добавляем на нее горизонтальную прямую \(y = 0\) (слой geom_hline):

ggplot(data = english, aes(x = WrittenFrequency, y = residuals)) +
  geom_point() + geom_hline(yintercept = 0, colour = 'red')

Что по такому графику мы можем определить? Во-первых, если какие-то закономерности в остатках. Если они есть, что свидетельствует о нарушении одного из условий модели, то точки «выстраиваются» на графике в определенном порядке, неслучайно, например, так:

Пример идеальной ситуации, когда остатки разбросаны случайным образом:

ansmod2 <- lm(data = anscombe, y1 ~ x1)
anscombe$res2 = ansmod2$residuals
ggplot(data = anscombe, aes(x = x1, y = res2)) +
  geom_point() + geom_hline(yintercept = 0, colour = 'red')

В нашем случае точек много, но, в целом, видно что при больших значениях WrittenFrequency точки менее разбросаны вокруг \(y=0\), поэтому нельзя говорить, что в распределении остатков нет никаких закономерностей. Второе условие не выполняется. По той же причине не выполняется и четвертое условие: разброс точек больше в левой части графика, чем в правой (посмотрите, как точки удалены от прямой \(y=0\)), поэтому дисперсию остатков нельзя считать постоянной. Третье условие более-менее выполняется: математическое ожидание остатков можно считать равным нулю, поскольку прямая \(y= 0\) лежит примерно посередине получившегося облака точек.

Для проверки пятого условия достаточно построить нормальную вероятностную бумагу для остатков (residuals) или применить один из критериев для проверки нормальности.

ggplot(data = english, aes(sample = residuals)) + stat_qq() + stat_qq_line()

На нормальное распределение не похоже. Проверим формально:

shapiro.test(english$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  english$residuals
## W = 0.96964, p-value < 2.2e-16

Гипотеза о нормальности распределения остатков отвергается.

Итоги: можем заключить, что качество нашей модели оставляет желать лучшего: \(R^2\) довольно низкий, условия Гаусса-Маркова не выполняются.