Практикум

Условие

Загрузите базу данных из исследования психологических детерминант волонтерства. Описание базы данных можно найти здесь.

vol <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Cowles.csv")

  1. Выведите описательные статистики для переменной экстраверсии. Есть ли в переменной нетипичные значения? Если есть, удалите их из анализа.

  2. Проверьте переменную экстраверсии на нормальность визуально: построив нормальную вероятностную бумагу. Предложите интерпретацию.

  3. Проверьте переменную экстраверсии на нормальность, используя критерии Колмогорова-Смирнова и Шапиро-Уилка. Какой вывод Вы сделаете относительно нулевой гипотезы?

  4. Выведите описательные статистики для переменной экстраверсии отдельно для женщин и для мужчин, используя пакет psych. Различаются ли средние значения по группам?

  5. Постройте гистограммы по группам для переменной экстраверсии отдельно для женщин и для мужчин.

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

  7. Проверьте гипотезу о равенстве среднего уровня экстраверсии среди мужчин и женщин. Каким критерием Вы будете пользоваться? Обоснуйте свой выбор. Какой вывод Вы можете сделать?

Решение

Сразу откроем пакеты, которые нам могут потребоваться в работе: psych для описательных статистик и ggplot2 для графиков.

library(psych)
library(ggplot2)
  1. Выведем описательные статистики для переменной экстраверсии.

Первый способ:

summary(vol$extraversion)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   10.00   13.00   12.37   15.00   23.00

Второй способ (с помощью пакета psych):

describe(vol$extraversion)
##    vars    n  mean   sd median trimmed  mad min max range  skew kurtosis
## X1    1 1421 12.37 3.89     13   12.43 4.45   2  23    21 -0.15    -0.34
##     se
## X1 0.1

Чтобы узнать, есть ли нетипичные значения, надо построить график ящика с усами. Судя по графике, такие значения в выборке присутствуют.

boxplot(vol$extraversion, ylab = "Extraversion", col = "salmon")

Удалим нетипичные значения, рассчитав границы типичных значений. Обратите внимание на то, что неравенство нестрогое.

\[[x_{0.25}-1.5\times IQR; x_{0.75}+1.5\times IQR]\]

lower_bound <- quantile(vol$extraversion, 0.25)-1.5*IQR(vol$extraversion)
upper_bound <- quantile(vol$extraversion, 0.75)+1.5*IQR(vol$extraversion)
vol_no_out <- subset(vol, extraversion >= lower_bound & extraversion <= upper_bound)

Нетипичные наблюдения успешно отфильтрованы из базы:

boxplot(vol_no_out$extraversion, ylab = "Extraversion", col = "salmon")

Дальше будем продолжать работать с очищенной базой.

  1. Построим нормальную вероятностную бумагу. График приблизительно похож на нормальный, не считая того, что в выборке много повторяющихся значений. Действительно, переменная экстраверсии принимает только целочисленные значения.
qqnorm(vol_no_out$extraversion)
qqline(vol_no_out$extraversion)

  1. Проверим на нормальность, используя статистические критерии. Оба теста свидетельствуют, что нулевую гипотезу о равенстве распределения нормальному нужно отвергнуть (p-value очень мало).
ks.test(vol_no_out$extraversion, "pnorm",
        mean = mean(vol_no_out$extraversion, na.rm = T),
             sd = sd(vol_no_out$extraversion, na.rm = T))
## Warning in ks.test(vol_no_out$extraversion, "pnorm", mean =
## mean(vol_no_out$extraversion, : ties should not be present for the
## Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  vol_no_out$extraversion
## D = 0.077021, p-value = 1.023e-07
## alternative hypothesis: two-sided
shapiro.test(vol_no_out$extraversion)
## 
##  Shapiro-Wilk normality test
## 
## data:  vol_no_out$extraversion
## W = 0.98882, p-value = 6.007e-09
  1. Средние значения экстраверсии у женщин и мужчин мало отличаются друг от друга и равны 12.45 и 12.35 соответственно.
describeBy(vol_no_out$extraversion, group = vol_no_out$sex)
## 
##  Descriptive statistics by group 
## group: female
##    vars   n  mean   sd median trimmed  mad min max range  skew kurtosis
## X1    1 776 12.45 3.76     13    12.5 4.45   3  22    19 -0.11    -0.37
##      se
## X1 0.13
## -------------------------------------------------------- 
## group: male
##    vars   n  mean   sd median trimmed  mad min max range  skew kurtosis
## X1    1 639 12.35 3.94     13    12.4 4.45   3  22    19 -0.14    -0.47
##      se
## X1 0.16
  1. Построим гистограммы отдельно для мужчин и для женщин.

Первый способ. Создадим отдельные датафреймы для подвыборок. Лайфхак: прежде чем фильтровать базу данных по категориальной переменной, посмотрите, как она закодирована! Если это факторная переменная (цифровое значение имеет текстовый лейбл), то фильтровать нужно по текстовому значению. Проверять можно, например, вот так:

table(vol_no_out$sex)
## 
## female   male 
##    776    639

Отлично, теперь мы видим, какие значения принимает переменная. А еще можно было просто открыть таблицу (командой View(vol_no_out) или нажав на значок таблицы в Global Environment) и посмотреть.

Фильтруем:

females <- subset(vol_no_out, sex == "female")
males <- subset(vol_no_out, sex == "male")

Строим гистограммы. Графики хоть и близки к нормальным распределениям, но все-таки видны расхождения (например, левые хвосты тяжелее, чем правые).

hist(females$extraversion, freq = F, col = "salmon", 
     main = "Extraversion among females", xlab = "Extraversion")
curve(dnorm(x, mean = mean(females$extraversion, na.rm = T), sd = sd(females$extraversion, na.rm = T)), add = T)

hist(males$extraversion, freq = F, col = "salmon",
     main = "Extraversion among males", xlab = "Extraversion")
curve(dnorm(x, mean = mean(males$extraversion, na.rm = T), sd = sd(males$extraversion, na.rm = T)), add = T)

Второй способ.

Нам потребуется пакет ggplot2. Обратите внимание, что аргументом binwidth я настраиваю длину столбца, равную 2 (баллам экстраверсии), чтобы график получился более информативным.

ggplot(data = vol_no_out, aes(x = extraversion, y = ..density..)) + 
  geom_histogram(fill = "coral", col = "black", binwidth = 2) + 
  facet_grid(~sex)

  1. Проверим на нормальность по подвыборкам.

Женщины:

ks.test(females$extraversion, "pnorm",
        mean = mean(females$extraversion, na.rm = T),
        sd = sd(females$extraversion, na.rm = T))
## Warning in ks.test(females$extraversion, "pnorm", mean =
## mean(females$extraversion, : ties should not be present for the Kolmogorov-
## Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  females$extraversion
## D = 0.078947, p-value = 0.0001259
## alternative hypothesis: two-sided
shapiro.test(females$extraversion)
## 
##  Shapiro-Wilk normality test
## 
## data:  females$extraversion
## W = 0.98913, p-value = 1.65e-05

Мужчины:

ks.test(males$extraversion, "pnorm",
        mean = mean(males$extraversion, na.rm = T),
        sd = sd(males$extraversion, na.rm = T))
## Warning in ks.test(males$extraversion, "pnorm", mean =
## mean(males$extraversion, : ties should not be present for the Kolmogorov-
## Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  males$extraversion
## D = 0.074478, p-value = 0.001668
## alternative hypothesis: two-sided
shapiro.test(males$extraversion)
## 
##  Shapiro-Wilk normality test
## 
## data:  males$extraversion
## W = 0.98727, p-value = 2.329e-05

Нормальности нет ни в одной группе. Гипотеза о равенстве нормальности отвергается согласно каждому из тестов (p-value очень мало).

  1. Так как наши данные не имеют нормального распределения, мы будем пользоваться непараметрическим тестом Вилкоксона. Он проверяет нулевую гипотезу о равенстве распределений в двух независимых группах.
wilcox.test(vol_no_out$extraversion ~ vol_no_out$sex)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  vol_no_out$extraversion by vol_no_out$sex
## W = 251160, p-value = 0.672
## alternative hypothesis: true location shift is not equal to 0

Нулевая гипотеза не отвергается (p-value = 0.672). Следовательно, мы можем сделать вывод о том, что распределения экстраверсии у мужчин и женщин одинаков.

Семинар

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

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

Корреляционный анализ всегда стоит начинать с построения диаграмм рассеяния, о которых мы говорили на занятии 7.

Давайте откроем базу нашего психологического эксперимента и визуализурем взаимосвязь между ростом и баллом ЕГЭ по математике.

psych <- read.csv2("http://math-info.hse.ru/f/2018-19/psych-ms/data/psych_survey.csv")
plot(psych$height, psych$math, xlab = "Рост", ylab = "Балл ЕГЭ по математике")
abline(lm(psych$math ~ psych$height))

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

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

Рассчет

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

Для расчета коэффициентов корреляции используется команда cor.test(x, y, method). По дефолту аргумент method принимает значение "person".

Внизу выдачи, после заглавия cor, приводится оценка коэффициента корреляции. Вверху выдачи приведены значения статистики критерия и p-value. Мы видим, что несмотря на то, что значение коэффициента корреляции Пирсона равно \(-0.148\), это значение статистически не отличается от нуля: p-value равно \(0.1293\). Значит, связи между ростом и баллом ЕГЭ по математике нет.

cor.test(psych$height, psych$math, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  psych$height and psych$math
## t = -1.529, df = 104, p-value = 0.1293
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.32970374  0.04371827
## sample estimates:
##        cor 
## -0.1482737

Выбросы

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

psych_no_out <- subset(psych, height > 150)
cor.test(psych_no_out$height, psych_no_out$math, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  psych_no_out$height and psych_no_out$math
## t = -0.94979, df = 102, p-value = 0.3445
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2811486  0.1007753
## sample estimates:
##         cor 
## -0.09363009

Теперь значение коэффициента корреляции стало меньше (\(-0.09\)), а p-value больше (\(0.324\)).

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

Рассчет

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

Давайте рассчитаем коэффициент корреляции Спирмена для роста и балла ЕГЭ по математике. Значение коэффициента очень мало, и нулевая гипотеза о равенстве коэффициента нулю не отвергается.

cor.test(psych$height, psych$math, method = "spearman")
## Warning in cor.test.default(psych$height, psych$math, method = "spearman"):
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  psych$height and psych$math
## S = 207050, p-value = 0.6604
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.04316801

Теперь переоценим его для исключенного выброса.

cor.test(psych_no_out$height, psych_no_out$math, method = "spearman")
## Warning in cor.test.default(psych_no_out$height, psych_no_out$math, method
## = "spearman"): Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  psych_no_out$height and psych_no_out$math
## S = 192590, p-value = 0.7827
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.02737379

Мы видим, что значение коэффициента почти не поменялось.

Нелинейная взаимосвязь

Давайте посмотрим на примере, как коэффициенты корреляции Пирсона и Спирмена реагируют на нелинейные взаимосвязи.

Сгенерируем выборку из стандартного нормального распределения объемом в 150 наблюдений.

sample1 <- rnorm(150, mean = 10)

Создадим вторую выборку, нелинейно преобразовав ее из первой.

sample2 <- sample1^16 + sample1^4

Посмотрим на взаимосвязь:

plot(sample1, sample2)

Рассчитаем коэффициент корреляции Пирсона:

cor.test(sample1, sample2, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sample1 and sample2
## t = 10.298, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5419351 0.7306878
## sample estimates:
##       cor 
## 0.6460834

Несмотря на то, что вторая выборка функционально зависит от первой, коэффициент корреляции Пирсона равен 0.64.

Зато коэффициент корреляции Спирмена равен 1:

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

Тренировочное задание

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

depression <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Ginzberg.csv")

  1. Выведите описательные статистики для переменных adjdep (шкала депрессивности) и adjsimp (шкала, измеряющая склонность респондента видеть жизнь в упрощенных черно-белых тонах). Какие значения принимают переменные? В каких шкалах они измерены?

  2. Нарисуйте гистограммы для каждой из переменных. Установите ширину столбца, равную 0.2.

Подсказка. Команда hist() позволяет выставить только количество столбцов. Чтобы настроить ширину столбца, нужно воспользоваться пакетом ggplot2. За ширину столбца отвечает параметр binwidth в слое geom_histogram()

  1. Нарисуйте диаграмму рассеяния между переменными adjdep и adjsimp (отложите adjdep по оси OY). Что Вы можете сказать о взаимосвязи между переменными?

  2. Рассчитайте коэффициент корреляции и проверьте его на значимость. Какими коэффициентами Вы можете воспользоваться?