Кто не ворчлив, тот весел?

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

Загрузим базу данных:

df <- read.csv("https://nvasilenok.github.io/hse/games/chip_and_dale.csv")

Выведем описательные статистики сразу для всех переменных в базе данных, чтобы предположить, есть ли в выборке потенциальные выбросы. На данном этапе мы будем обращать внимание на минимумы и максимумы. Обратите внимание, что для качественных признаков R не вычисляет описательные статистики, а выводит распределение респондентов по категориям ответов.

summary(df)
##     group       gender      specialization resc_ranger_ch resc_ranger_now
##  191   :16   female:42   analysis  :29     Chip  :17      Chip  :23      
##  192   :17   male  :32   management:45     Dale  :12      Dale  :15      
##  193   :15                                 Gadget:31      Gadget:18      
##  194   :14                                 Monty : 6      Monty : 7      
##  195   :10                                 Other : 4      Other : 6      
##  prepod: 2                                 Zipper: 4      Zipper: 5      
##       dale             chip             monty           gadget     
##  Min.   :  0.00   Min.   :   5.00   Min.   :  0.0   Min.   :  0.0  
##  1st Qu.: 42.00   1st Qu.:  50.00   1st Qu.: 50.0   1st Qu.: 20.0  
##  Median : 67.50   Median :  63.50   Median : 72.5   Median : 55.0  
##  Mean   : 61.19   Mean   :  79.34   Mean   : 69.3   Mean   : 49.7  
##  3rd Qu.: 80.00   3rd Qu.:  90.00   3rd Qu.:100.0   3rd Qu.: 70.0  
##  Max.   :100.00   Max.   :1000.00   Max.   :101.0   Max.   :100.0  
##      zipper      
##  Min.   :  0.00  
##  1st Qu.: 42.50  
##  Median : 70.00  
##  Mean   : 59.57  
##  3rd Qu.: 78.92  
##  Max.   :100.00

Неестественно большое значение 1000 есть только в переменной chip, которую мы рассматривали на прошлом семинаре, во всех остальных переменных ничего подобного нет. Давайте исключим это наблюдение из выборки, сохранив новую выбору в новую таблицу данных:

df2 <- subset(df, chip < 1000)

Для рассчета коэффициентов корреляции и проверки их на значимость используется функция cor.test(). В качестве аргументов она принимает две переменные (они должны содержать одинаковое количество элементов), а также метод вычисления корреляции: method = "pearson" или method = "spearman". Если Вы не указываете аргумент method, то по дефолту рассчитывается коэффициент корреляции Пирсона.

Сначала оценим коэффициент корреляции Пирсона на первоначальной базе с выборосом. Внизу выдача содержит выборочную оценку коэффициента cor; вверху приведены результаты проверки гипотезы \(H_0:~Cor(X,Y)=0\) против двухсторонней альтернативы. Вычислены наблюдаемое значение статистики критерия t, количество степеней свободы df (вспомним, что для коэффициента корреляции Пирсона \(df=n-2\)) и p-value.

cor.test(df$dale, df$chip)
## 
##  Pearson's product-moment correlation
## 
## data:  df$dale and df$chip
## t = -2.8102, df = 72, p-value = 0.006372
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5065052 -0.0925424
## sample estimates:
##       cor 
## -0.314393

Мы видим, что \(R=-0.31\), что говорит о слабой отрицательной взаимосвязи, а \(\text{p-value}=0.006\), что говорит об отвержении гипотезы о незначимости. Следовательно, мы могли бы сделать вывод, что веселость и ворчливость связаны слабо, но отрицательно среди студентов 1 курса ОП “Политология”.

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

cor.test(df2$dale, df2$chip)
## 
##  Pearson's product-moment correlation
## 
## data:  df2$dale and df2$chip
## t = -2.1044, df = 71, p-value = 0.03888
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.44743316 -0.01296459
## sample estimates:
##        cor 
## -0.2423091

Действительно, новый коэффциент корреляции стал слабее \(R=-0.24\). Более того, \(\text{p-value}\) в разы уменьшился! Тем не менее, даже с новым \(\text{p-value}=0.04\) мы отвергаем нулевую гипотезу. Следовательно, студенты 1 курса Политологии все-таки чем менее ворчливы, тем более веселы.

Теперь посмотрим, как ведет себя коэффициент корреляции Спирмена. Оценим его на первоначальной базе:

cor.test(df$dale, df$chip, method = "spearman")
## Warning in cor.test.default(df$dale, df$chip, method = "spearman"): Cannot
## compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  df$dale and df$chip
## S = 83837, p-value = 0.03813
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.2415694

Структура выдачи похожа на коэффициент корреляции Пирсона: мы имеем выборочную оценку rho, а также статистику критерия S (обратите внимание, что R выводит сумму квадратов расстояний рангов) и p-value.

cor.test(df2$dale, df2$chip)
## 
##  Pearson's product-moment correlation
## 
## data:  df2$dale and df2$chip
## t = -2.1044, df = 71, p-value = 0.03888
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.44743316 -0.01296459
## sample estimates:
##        cor 
## -0.2423091

Мы видим, что \(\hat{\rho}=-0.24\), а \(\text{p-value}=0.04\). Следовательно, мы отвергаем нулевую гипотезу об отсутствии связи на 5% уровне значимости: связь значимая и отрицательная, но слабая. Результаты оценивания коэффициента корреляции Спирмена на таблице с выбросом очень похожи на коэффициент корреляции Пирсона без выброса!

(Не обращайте внимание на warning: он предупреждает, что в выборке есть одинаковые значения, и что R приходится работать с одинаковыми рангами.)

А как же сыр?

Теперь давайте узнаем, есть ли связь между веселостью dale и любовью к сыру monty. Первый шаг при вычислении корреляции всегда состоит в визуализации взаимосвязи: она позволяет предположить, какой будет связь, а также обнаружить потенциальные выбросы. Построим диаграмму рассеяния:

plot(df$dale, df$monty)

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

cor.test(df$dale, df$monty)
## 
##  Pearson's product-moment correlation
## 
## data:  df$dale and df$monty
## t = 0.73974, df = 72, p-value = 0.4619
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1445174  0.3092117
## sample estimates:
##        cor 
## 0.08684929

Действительно: мы видим, что оценка коэффициента корреляции оказалась очень близка к нулю \(R=0.09\), а очень высокий \(\text{p-value}=0.4619\) говорит о невозможности отвергнуть нулевую гипотезу. Судя по всему, связи между веселостью и любовью к сыру среди студентов 1 курса ОП “Политология” нет.

Линейная и монотонная взаимосвязи

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

Сгенерируем выборку из 1000 наблюдений и сохраним в вектор s, а затем преобразуем ее двумя способами. Сначала преобразуем линейно, сохранив в вектор s2. Затем преобразуем ее нелинейно с помощью мотононной функции \(f(x)=x^3\) и сохраним в вектор s3. Обратите внимание, что в обоих случаях трансфорации исходной выборки не предполагают никакой неопределенности и вариативности: взаимосвязь функциональна. И уже на этом этапе мы можем предсказать, что коэффициент корреляции Пирсона между s и s2 будет равен 1, и коэффициент корреляции Спирмена между s и s3 также будет равен 1. Вопрос в том, как будет вести себя коэффициент корреляции Писона при нелинейной взаимосвязи.

s <- rnorm(1000)
s2 <- s*10+10
s3 <- s^3

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

plot(s, s2)

plot(s, s3)

Теперь рассчитаем коэффициенты корреляции для каждой пары векторов. Оба коэффициента корреляции сумели поймать функциональную зависимость между s и s2: оба коэффициента корреляции равны 1.

cor.test(s, s2, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  s and s2
## t = 2.12e+09, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  1 1
## sample estimates:
## cor 
##   1
cor.test(s, s2, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  s and s2
## S = 0, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho 
##   1

Как мы и ожидали, коэффициент корреляции Спирмена между s и s3 равен 1, Но коэффициент корреляции Пирсона не может выявить нелинейную функциональную завивимость: его значение достаточно велико, но, тем не менее, отлично от 1.

cor.test(s, s3, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  s and s3
## t = 35.515, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7184614 0.7733474
## sample estimates:
##       cor 
## 0.7471755
cor.test(s, s3, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  s and s3
## S = 0, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho 
##   1