21 сентября 2012

Классические методы статистики: коэффициент корреляции



Коэффициент корреляции (\(r\)) - очень удобный показатель степени взаимосвязи между двумя переменными. Он представляет собой безразмерную величину, которая изменяется от \(-1\) до \(+1\). При независимом варьировании переменных, когда связь между ними отсутствует, \(r = 0\). Чем сильнее связь, тем больше величина коэффициента корреляции. При этом положительные значения \(r\) указывают на положительную (= прямую) связь (т.е. при увеличении значений одной переменной в среднем возрастают значения и другой переменной), а отрицательные - на отрицательную (= обратную) связь (при возрастании одной переменной другая уменьшается).

Вычисление коэффициента корреляции в R я продемонстрирую на примере данных из своей статьи (Mastitsky 2012), которые можно напрямую загрузить в R с сайта figshare.

dat <- read.delim("http://figshare.com/media/download/98923/97987")
 
head(dat)
  Month     Lake Site ZMlength CAnumber
1   May Batorino   S3     14.9       36
2   May Batorino   S3     14.0       30
3   May Batorino   S3     13.0      331
4   May Batorino   S3     14.0      110
5   May Batorino   S3     12.0        4
6   May Batorino   S3     14.0      171

Источник: Mastitsky et al. 2008
Таблица dat содержит данные по уровню зараженности двустворчатого моллюска Dreissena polymorpha инфузорией-комменсалом Conchophthirus acuminatus (см. фото слева) в трех озерах Беларуси, различающихся по уровню трофности. Нас интересуют, в частности, две переменные: длина раковины моллюска (ZMlength, мм) и число обнаруженных в моллюлюске инфузорий (CAnumber). На приведенном ниже рисунке прослеживается положительная связь между этими двумя переменными (замечание: в рассматриваемом примере анализируются все данные, без разделения по озерам). Вопрос, однако, состоит в том, насколько сильна эта связь. Оценить ее поможет коэффициент корреляции.

Связь между длиной раковины дрейссены и интенсивностью инвазии Conchophthirus acuminatus

Начнем с коэффициента корреляции Пирсона (англ. Pearson correlation coefficient), который, как известно, рассчитывается по формуле:

\[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}} \]

Конечно, выполнять вычисления вручную нет необходимости. В R коэффициент корреляции Пирсона, равно как и другие коэффициенты, можно легко рассчитать при помощи функциий cor() и cor.test(). Различие между этими двумя функциями заключается в том, что cor() позволяет вычислить только сам коэффициент корреляции, тогда как cor.test() выполняет еще и оценку статистической значимости коэффициента, проверяя нулевую гипотезу о равенстве его нулю. Я предпочитаю использовать именно вторую функцию:


attach(dat)
 
cor.test(CAnumber, ZMlength)
 
 Pearson's product-moment correlation
 
data:  CAnumber and ZMlength 
t = 11.4964, df = 474, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.3935877 0.5343949 
sample estimates:
     cor 
0.466946 

Как видим, рассчитанный коэффициент корреляции Пирсона оказался равен 0.467. Несмотря на то, что он не очень высок, этот коэффициент статистически значимо отличается от нуля (p-value < 2.2e-16). Для нашего удобства, программа также автоматически вычислила 95%-ный доверительный интервал для полученного коэффициента корреляции (95 percent confidence interval: 0.394  0.534).


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

  • Обе анализируемые переменные распределены нормально
  • Связь между этими переменными линейна

Приведенный ниже рисунок показывает, что как минимум в отношении значений интенсивности инвазии условие нормальности распределения не выполняется:
 

Для исправления ситуации можно попробовать логарифмировать обе переменные, т.е. и ZMlength, и CAnumber:




Для преобразованных переменных коэффициент корреляции Пирсона составит:

cor.test(log(CAnumber+1), log(ZMlength))
 
 Pearson's product-moment correlation
 
data:  log(CAnumber + 1) and log(ZMlength) 
t = 21.5166, df = 474, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.6543961 0.7456953 
sample estimates:
      cor 
0.7029297

Неудивительно, что новое значение коэффициента корреляции значительно выросло (0.703 против 0.467). И все же, несмотря на логарифирование, значения интенсинвости инвазии C. acuminatus не подчиняются нормальному распределению. Показать это можно как графически, таки и при помощи, например, теста Шапиро-Уилка (подробнее о проверке данных на нормальность см. здесь):
 
shapiro.test(log(CAnumber+1))
 
 Shapiro-Wilk normality test
 
data:  log(CAnumber + 1) 
W = 0.9508, p-value = 1.734e-11


Для ненормально распределенных переменных, а также при наличии нелинейной связи между переменными, следует использовать непараметрический коэффициент корреляции Спирмена (англ. Spearman correlation coefficient). В отличие от коэффициента Пирсона, этот вариант коэффициента корреляции работает не с исходными значениями переменных, а с их рангами (формула при этом используется та же, что и для коэффициента Пирсона - см. выше). Для вычисления коэффициента Спирмена в R при вызове функции cor.test() необходимо воспользовать аргументом method со значением "spearman":

cor.test(CAnumber, ZMlength, method = "spearman")
 
 Spearman's rank correlation rho
 
data:  CAnumber and ZMlength 
S = 6574110, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0 
sample estimates:
      rho 
0.6342627 
 
Warning message:
In cor.test.default(CAnumber, ZMlength, method = "spearman") :
  Cannot compute exact p-values with ties
 





Коэффициент корреляции Спирмена составил 0.634 и оказался статистически значимым (Р << 0.001). Поскольку в данных имеют место значения с одинаковыми рангами ("связанные ранги", англ. ties), программа не смогла рассчитать точное Р-значение, о чем предупредила в сообщении "Warning message: ... Cannot compute exact p-values with ties". В связи с тем, что коэффициент корреляции Спирмена работает с рангами, любое преобразование исходных данных никак не сказывается на его значении. Например, после логарифмирования мы получим результат, идентичный предыдущему:

cor.test(log(CAnumber+1), log(ZMlength), method = "spearman")
 
 Spearman's rank correlation rho
 
data:  log(CAnumber + 1) and log(ZMlength) 
S = 6574110, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0 
sample estimates:
      rho 
0.6342627

Наконец, третий вариант коэффициента корреляции, который можно рассчитать при помощи функции cor.test(), имеет название коэффициент ранговой корреляции Кендалла, \( \tau \) (англ. Kendall's tau). Работает он следующим образом. Предположим, что у нас есть набор из парных наблюдений для двух переменных: \( (x_1, y_1), (x_2, y_2)  ... (x_n, y_n) \). Говорят, что две пары наблюдений \( (x_i, y_i) \) и \( (x_j, y_j) \) являются конкордантными, если имеется согласие между рангами соответствующих элементов этих пар, т.е. если
  • \( x_i > x_j \) и \( y_i > y_j \)
  • или  \( x_i < x_j \) и \( y_i < y_j \)
Коэффициент Кендалла рассчитывается следующим образом:

\[ \tau = \frac{n_{conc. pairs} - n_{discord. pairs}}{0.5n(n-1)} \]

где \( n_{conc. pairs}\) - число конкордантных пар, а \( n_{discord. pairs}\) - число дискордантных пар. Коэффициент Кендалла часто используют при анализе того, насколько хорошо согласуются результаты измерений, получаемые при помощи разных приборов, или результаты голосований экспертов по одному и тому же вопросу, и т.п.

В R коэффициент Кендалла можно вычислить так:

cor.test(CAnumber, ZMlength, method = "kendall")
 
 Kendalls rank correlation tau
 
data:  CAnumber and ZMlength 
z = 14.9307, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0 
sample estimates:
      tau 
0.4655232

Подробнее об аргументах функции cor.test() можно узнать из ее справочного файла, доступного по команде ?cor.test


12 комментариев :

Александр Марцинкевич комментирует...

Подскажите, а как можно сравнить rho Спирмена? Необходимо сделать заключение о различии/схожести двух коэффициентов корреляции. Для r Пирсона решение есть - через критерий Фишера, а для Спирмена?

Анонимный комментирует...

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

У меня возникли два небольших вопроса (не судите строго ламера).

1. При расчете различных коэффициентов и p-значений, программа R часто выдает ответы типа (на примере ваших результатов):
- p-value < 2.2e-16;
- p-value = 1.734e-11;
- p-value < 2.2e-16.
Подскажите пожалуйста, как понимать эти показатели статистической значимости? Что означает буква "е" и стоящий после неё "-" в контексте числовых обозначений?

2. Хотелось бы больше узнать о процедуре преобразования данных: в чем её суть, какие ещё способы преобразования существуют, кроме нахождения логарифмов?

3. При расчете коэффициента корреляции Пирсона -
cor.test(log(CAnumber+1), log(ZMlength)),
к одной из логарифмированых переменных (CAnumber) Вы добавляете единицу, а к другой переменной (ZMlength) - нет.
Для чего это делается?


Sergey Mastitsky комментирует...

1. Буква "е" в выражениях вроде "p-value < 2.2e-16" буквально означает "умножить на 10 в степени...". То есть 2.2e-16 = 2 * 10^(-16). Как Вы можете догадаться, получится очень и очень малое число, что указывает на очень высокую значимость статистического теста (в случае с коэф. корреляции мы проверяем гипотезу о равенстве его нулю; соответственно при таких низких р-значениях мы можем смело эту гипотезу отвергать).

2. Достаточно подробно об этом можно почитать в статье Википедии (англ. яз.):
https://en.wikipedia.org/wiki/Data_transformation_%28statistics%29
Об одном из примеров, где может потребоваться преобразование данных, (трансформация) я писал здесь: http://r-analytics.blogspot.de/2013/05/blog-post.html#.UeL7tMgmOHs
Крайне рекомендую также классическую статью: http://pegasus.cc.ucf.edu/~lni/sta6236/BoxCox1964.pdf

3. Единицу была добавлена к переменной, которая содержит нулевые значения. Причина проста - логарифм нуля не определен. Поскольку единица добавляется ко всем значением переменной, это совершенно никак не сказывается на конечном результате - величине коэффициента корреляции. Но при этом мы решаем проблему с log(0). Это распространенный подход для таких ситуаций. Часто можно встретить случаи, когда добавляют не 1, а какие-то другие небольшие величины, например, 0.001 или 0.5. Что именно добавлять, будет зависеть от свойст самой переменной. В примере, рассмотренном, в сообщении, 1 - логичный выбор, поскольку это наименьшее ненулевое значения для CAnumber.

Анонимный комментирует...

Сергей, благодарю Вас за оперативный и содержательный ответ, и вообще за этот проект - Вы делаете хорошее дело.
Спасибо, также, за источники - буду изучать.

Анонимный комментирует...

Подскажите пожалуйста, во входящих значениях должна отсутствовать автокорреляция при проведении теста на корреляцию между ними или ее наличие не важно? заранее спасибо.

Sergey Mastitsky комментирует...

Очень мудрено заданный вопрос, на который, соответственно, сложно ответить. Что Вы подразумеваете под "входящими значениями"? Вы хотите выяснить корреляцию между двумя векторами числовых значений или обнаружить автокорреляцию в одном конкретном временном ряду?

Анонимный комментирует...

Сергей, имеются два ряда (числовые векторы), в одном из них присутствует автокорреляция. Можно ли проводить тест на наличие между этими двумя рядами корреляции. Где-то попалась заметка (не могу найти), что коэфф. корреляции правильно показывает тесноту связи если в рядах отсутствует автокорреляция. Больше ни где не видела, вот и хотела узнать насколько это важно, отсутствие автокорреляции в рядах при поиске связи между ними?

Sergey Mastitsky комментирует...

Проводить тест на наличие корреляции, конечно, никто не запрещает. Но Вы правы в том, что обычный (t-тест) тест на значимость коэффициента корреляции (Пирсона) может выдать некорректное Р-значение. Вкратце, проблема в том, что при наличии выраженной автокорреляции стандартная ошибка соответствующей переменной будет оценена неверно (обычно завышена), что, в свою очередь, повлияет на t-статистику и соответствующее Р-значения. Решения могут быть разными, но все они будут сводиться к внесению поправки в оценки стандартных ошибок.
Обсуждение похожих проблем Вы можете найти на http://stats.stackexchange.com/
и в некоторых статьях, например: http://www.sandyliebhold.com/pubs/florence.PDF

Анонимный комментирует...

Ок! Спасибо Сергей большое.

Анонимный комментирует...

скажите пожалуйста,а если нужно найти корреляцию трех ранжировок?или более?как тогда?

Sergey Mastitsky комментирует...

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

Анонимный комментирует...

Тест Пирсона всегда дает меньший результат, чем тест Спирмена. Скажите, пожалуйста, с чем это связано? Можно ли на этом основании делать какой-либо вывод?

Отправить комментарий