Лекция №3: Корреляции

Георгий Мороз, Ольга Ляшевская, Илья Щуров

Курс «Анализ лингвистических данных: квантитативные методы и визуализация».

Школа лингвистики НИУ ВШЭ, магистерская программа «Теория языка и компьютерная лингвистика», 2015-16 учебный год.

Школа лингвистики | Все материалы курса | Исходные коды на Github

Данный материал доступен под лицензией CC BY-SA 4.0. При использовании обязательно упоминание авторов курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса.


Что такое корреляция

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

Рассмотрим датафрейм, в котором есть две переменные (x и y) и сколько-то наблюдений.

correlatedValue = function(x, r, mean=0) {
  r2 = r**2
  ve = 1-r2
  SD = sqrt(ve)
  e  = rnorm(length(x), mean=mean, sd=SD)
  y  = r*x + e
  return(y)
}
# the function above is from here: http://stats.stackexchange.com/a/38867

library(knitr)
N <- 200
K <- 5
x <- rnorm(N) + 2
y <- correlatedValue(x, 0.8, -1)
df <- data.frame(x,y)
rownames(df) <- 1:N
kable(df[1:K, ], row.names = T, caption = sprintf("В таблице приведены данные только по первым %i наблюдениям, но на самом деле их больше.", K))
В таблице приведены данные только по первым 5 наблюдениям, но на самом деле их больше.
x y
1 0.6149294 0.4720610
2 2.0383232 0.4046687
3 1.2369698 -0.1101671
4 2.2123061 0.7505535
5 3.4255380 2.0483225
library(ggplot2)

Как нам что-то узнать про связь между x и y? Начнём с картинки.

ggplot(df)+geom_point(aes(x,y))+ggtitle('Рис. 1')

Здесь изображён scatter plot (точечная диаграмма или диаграмма рассеяния) для нашего датафрейма. Каждой точке соответствует одно наблюдение: горизонтальная координата равна значению переменной x, а вертикальная — переменной y.

Из картинки видно, что облако точек вытянуто и наклонено. Это значит, что большим значениям переменной x соответствуют большие значения y. Это не строгое правило: бывают такие пары точек, что координата x больше у первой, а y у второй, но в среднем это верно. Например, если выбрать все точки с x-координатой от 3 до 4, то среднее значение y-координат для них будет больше, чем аналогичное среднее значение для точек с x-координатой от 0 до 1.

Посмотрим ещё на несколько картинок:

plot_cor_example = function(x, r, main, size = 2) {
  y <- correlatedValue(x, r)
  df <- data.frame(x,y)
  list(ggplot(df)+geom_point(aes(x, y), size=size)+ggtitle(main), df)
}

library(cowplot)

x <- rnorm(N)+2
A <- plot_cor_example(x, 0.6, main = 'Рис. 2')
B <- plot_cor_example(x, 0.95, main = 'Рис. 3')
C <- plot_cor_example(x, 0, main = 'Рис. 4')
D <- plot_cor_example(x, -0.7, main = 'Рис. 5')
plot_grid(A[[1]],B[[1]],C[[1]],D[[1]])


Как можно их охарактеризовать? На рисунке 2 зависимость проявляется гораздо слабее, по сравнению с рисунком 1, но её по-прежнему можно увидеть. На рисунке 3 зависимость проявляется очень явно: точки выстроились практически вдоль прямой с ненулевым наклоном. На рисунке 4 никакой ярко выраженной зависимости не видно: облако точек не имеет никакого явного наклона. На рисунке 5 зависимость есть, но она обратная: чем больше x, тем меньше y.

Картинки — это хорошо, но нам бы хотелось чего-то более осязаемого: делать выводы по одной только картинке опасно — вдруг вы видите в ней то, чего не увидят другие? Поэтому нам хочется придумать какую-то числовую характеристику, которая бы позволила отвечать на вопрос о наличии и характере зависимости. Такая величина называется коэффициентом корреляции или иногда просто корреляцией. Коэффициенты корреляции бывают разные. Мы начнём с корреляции Пирсона.

Корреляция Пирсона

Вернёмся к картинке на рисунке 1. Для начала сдвинем систему координат в точку \( (\bar x, \bar y) \), где \( \bar{x}\) и \( \bar y\)— средние значения соответствующих переменных. Иными словами, будем вместо точек с координатами (xi, yi) рассматривать точки с координатами \( (x_i-\bar x, y_i-\bar y) \). Получится такая картинка:

pearsonize <- function(df, legend = T, size=3) {
  df$new_x <- df$x - mean(df$x)
  df$new_y <- df$y - mean(df$y)
  df$sign <- as.factor((df$new_x * df$new_y)>0)
  if (legend) {
    legend = scale_color_discrete(name='new_x*new_y', breaks=c(T,F), labels=c(">0","<0"))
  }
  else{
    legend = scale_color_discrete(guide=F)
  }
  ggplot(df) +
    geom_point(aes(new_x,new_y, colour=sign), size=size) +
    legend + 
    geom_vline(xintercept=0, colour='gray', linetype='longdash') + 
    geom_hline(yintercept=0, linetype='longdash', colour='gray')
}
pearsonize(df) + ggtitle("Рис. 6")

Для точки с номером i вычислим теперь произведение \( (x_i-\bar x) (y_i-\bar y) \). Заметим, что для тех точек, которые лежат в правом верхнем и левом нижнем квадрантах это произведение положительно, а для тех, которые лежат в левом верхнем и правом нижнем квадрантах оно отрицательное. Наш график вытянут в направлении «из левого нижнего угла в правый верхний», поэтому точек с положительным произведением больше, чем с отрицательным. Если теперь сложить все такие произведения, то получающееся число будет положительным. Если поделить его на число точек n, получится величина, называемая выборочной ковариацией:


$$ \mathop{\mathrm{cov}} (x,y)=\frac{1}{n}\sum_{i=1}^{n} (x_i-\bar x)(y_i-\bar y)=\frac{(x_1 - \bar x)(y_1-\bar x)+(x_2-\bar x)(y_2 - \bar y)+\ldots+(x_n-\bar x)(y_n - \bar y)}{n}, $$

Посмотрим теперь на картинки с рисунков 2—5.

plot_grid(
  pearsonize(A[[2]], legend=F, size=2)+ggtitle("Рис. 2'"),
  pearsonize(B[[2]], legend=F, size=2)+ggtitle("Рис. 3'"),
  pearsonize(C[[2]], legend=F, size=2)+ggtitle("Рис. 4'"),
  pearsonize(D[[2]], legend=F, size=2)+ggtitle("Рис. 5'")
)

На рисунке 2' бирюзовые точки с положительным произведением по-прежнему преобладают над красными точками с отрицательным произведением, но последних по сравнению с рисунком 6 гораздо больше. На рисунке 3' наоборот красных точек почти нет. На рисунке 4' красных и бирюзовых точек почти поровну, а на рисунке 5' красных точек больше, чем бирюзовых. Если мы посчитаем выборочную ковариацию для каждой из этих картинок, было бы разумно ожидать, что для рисунка 3' она окажется самой большой, для рисунка 4' окажется близкой к нулю, а для рисунка 5' вообще отрицательной.

Кажется, что ковариация — это то, что нам нужно. Но не будем торопиться: у неё есть одна большая проблема. Представим себе, что мы изменили единицу измерения величины x: например, вместо метров стали использовать миллиметры. В этом случае все xi увеличились в тысячу раз! Если посмотреть на формулу для ковариации, видно, что в этом случае и ковариация увеличится в тысячу раз. Но ведь на характере связи изменение единицы измерения не должно отражаться — связь между переменными не исчезнет, не усилится и не ослабнет. Значит, ковариация — это не совсем то, что нам нужно.

Чтобы исправить этот недостаток ковариации, нужно разделить её на произведение стандартных отклонений переменных x и y. Получится такая штука:


$$\mathop{\mathrm{cor}}(x,y)=\frac{\mathop{\mathrm{cov}}(x,y)}{SD(x)\cdot SD(y)},$$

где SD — это стандартное отклонение. Теперь если все значения переменной x увеличить в 1000 раз, то в 1000 раз увеличится \( \mathop{\mathrm{cov}}(x,y) \) (посмотрите внимательно на формулу с определением ковариации и проверьте, что это так) и одновременно в 1000 раз увеличится стандартное отклонение SD(x), а значит дробь не изменится.

Получающаяся величина называется коэффициентом корреляции Пирсона. Этот коэффициент всегда находится на отрезке [ − 1, 1]. Нетрудно показать (если вы любите возиться с формулами, сделайте это как упражнение), что коэффициент корреляции равен 1 в том и только в том случае, когда точки на картинке располагаются ровно вдоль прямой, причём эта прямая имеет положительный наклон (если задать её функцией y = kx + b, то k > 0). А если прямая имеет отрицательный наклон (k < 0), то коэффициент корреляции равен −1. Чем ближе точки к прямой, тем больше по модулю коэффициент корреляции.

Построим ещё несколько картинок и посчитаем для них коэффициент корреляции Пирсона:

corrs <- list(-1, -0.9, -0.7, -0.4, 0, 0.4, 0.7, 0.9, 1)
x <- rnorm(N)+2
do.call(plot_grid, lapply(corrs, function(r) {plot_cor_example(x,r,main=as.character(r), size=1)[[1]]}))

Внимание! Коэффициент корреляции Пирсона реагирует только на близость точек к прямой, но не на её наклон. Это следует ровно из того факта, что он не меняется при умножении всех значений x или всех значений y на какое-то число. Единственное, что можно сказать о наклоне прямой, зная коэффициент корреляции — это положительной он или отрицательный. Для того, чтобы оценивать наклон, используются другие инструменты — регрессии. О них мы поговорим в следующий раз.

Когда Пирсон не помогает: нелинейная зависимость

Сравним два рисунка.

N <- 50
x <- abs(rnorm(N))
y <- sin(x**(1/10)*1.3)
df <- data.frame(x,y)
A <- ggplot(df) + geom_point(aes(x,y)) + ggtitle("Рис. 7a")
B <- plot_cor_example(rnorm(N), cor(x,y), "Рис. 7б")
plot_grid(A,B[[1]])

Совершенно понятно, что связь между переменными x и y гораздо более ярко выражена для данных, изображённых на рисунке 7a. Однако, коэффициент корреляции Пирсона для этих двух рисунков одинаковый (и равен примерно 0.79). Почему так происходит? Дело в том, что этот коэффициент меряет близость точек к некоторой прямой и измеряет уровень линейной связи между двумя переменными. Если эта связь не является линейной, то есть выражается не формулой вида y = kx + b, а какой-то другой, то точки располагаются вдоль какой-то другой линии и коэффициент корреляции Пирсона уменьшается.

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

Совсем другое дело! После логарифмирования коэффициент корреляции равен 0.99.

Корреляция для порядковых шкал

Можно подойти к рисунку 7a с другой стороны. Здесь, очевидно, выполняется следующее правило: для любых двух точек, та, которая находится правее, находится также и выше. Иными словами, если x-координата больше, то и y-координата больше. Обратите внимание: чтобы проверить это условие, нам вообще не нужно ничего знать про координаты, кроме того, в каком порядке они идут друг относительно друга. Иными словами, мы можем перейти от абсолютной шкалы к порядковой!

Для порядковых шкал есть два распространённых способа считать корреляцию: корреляция Спирмена и корреляция Кенделла. Рассмотрим их по порядку.

Корреляция Спирмена

s <- c(7, 5, 9, 3, 1, 0)

Рассмотрим какую-нибудь выборку, то есть просто ряд чисел. Например, 7, 5, 9, 3, 1, 0. Можно упорядочить эти числа по возрастанию: 0, 1, 3, 5, 7, 9. Теперь можно каждому числу присвоить ранг — это просто его номер после упорядочивания:

x rank
0 1
1 2
3 3
5 4
7 5
9 6

Заменим теперь в исходной выборке каждое число на его ранг. Получится новый ряд 5, 4, 6, 3, 2, 1 — он называется «вариационным рядом» для исходной выборки. Если в исходном ряду есть несколько одинаковых чисел, то они получают одинаковый ранг, равный среднему арифметическому их номеров при сортировке. Например, вариационный ряд для выборки 1 2 2 3 равен 1 2.5 2.5 4. Числа в вариационном ряду идут в том же порядке, как в исходной выборке, но вся остальная информация об элементах выборки полностью утеряна. Вот и хорошо. Этого мы и добивались.

Пусть теперь даны две выборки x и y. Вычислим от каждой из них вариационный ряд. Получим два набора целых чисел, они обозначаются \( \mathop{\mathrm{rank}}(x) \) и \( \mathop{\mathrm{rank}}(y) \). Посчитаем от них коэффициент корреляции Пирсона. То, что получится, называется коэффициентом корреляции Спирмена от исходной пары выборок x, y. Он обычно обозначается буквой \( \rho \).

Посмотрим ещё раз на рисунок 7a. Занумеруем для простоты точки в порядке возрастания координаты x. Тогда вариационный ряд для x будет самым простым: 1, 2, 3, .... А каким будет вариационный ряд для y? Точно таким же! Значит, коэффициент корреляции Пирсона для этой картинки равен 1. Звучит вполне разумно.

Корреляция Кенделла

Эта корреляция в некотором смысле даже проще, чем корреляция Пирсона. Давайте возьмём пару точек A и B: пусть точка A имеет координаты (xA, yA), а точка B координаты (xB, yB). Будем считать, что мы обозначили через A ту точку, у которой x-координата меньше, чем у другой. Тогда возможно два варианта: либо yA > yB, либо yA < yB. (Конечно, возможен ещё вариант, что координаты совпадают, но такое встречается редко и мы не будем углубляться в детали.) Скажем, что в первом случае пара точек \( \{A, B\} \) упорядочена правильно, а во втором случае неправильно.

Коэффициент корреляции Кенделла вычисляется так: возьмём все возможные пары точек на нашей картинке (то есть все возможные пары наблюдений). Всего их будет N штук. (Из комбинаторики мгновенно следует, что N = n(n − 1)/2, где n — число наблюдений, но это сейчас не очень важно.) Пусть среди них K пар упорядочено правильно и L пар неправильно. (В сумме K + L = N.) Тогда коэффициентом корреляции Кенделла (обозначается обычно \( \tau \) называется отношение


$$\tau = \frac{K-L}{N}$$

Вернёмся к рисунку 7a. Нетрудно видеть, что на нём любая пара точек упорядочена правильно. Значит, в формуле для \( \tau \) число K = N и L = 0, а следовательно \( \tau=1 \). Тоже согласуется с реальностью.

Как вычислить корреляцию в R

Очень просто. Для этого используется функция cor.

x <- c(8, 1, 2, 9, 3)
y <- c(12, 1, 2, 8, 4)
cor(x,y)
## [1] 0.9138735

По умолчанию она вычисляет корреляцию Пирсона, но можно попросить вычислить и Кенделла, и Спирмена:

cor(x,y, method='k') # Kendall
## [1] 0.8
cor(x,y, method='s') # Spearman
## [1] 0.9