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

Коэффициент корреляции (\(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