Семинар

В конспект этого семинара входят три важные для нашего курса темы:

  1. Визуальная проверка на нормальность: нормальная вероятностная бумага (Q-Q plot).
  2. Статистическая оценка как случайная величина.
  3. Проверка статистических гипотез: равенство распределения нормальному и равенство средних.

Продолжаем работать с базой данных предыдущих занятий psych_survey.csv

Визуальная проверка на нормальность: нормальная вероятностная бумага

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

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

hist(df$height, breaks = 20, freq = F,
     col = "navajowhite", xlab = "Height", main = "Histogram of height")
curve(dnorm(x, mean(df$height), sd(df$height)), 
      add = T)

Есть другой способ графически сравнивать распределение с нормальным. Для этого строится график нормальной вероятностной бумаги (q-q plot, quantile-quantile plot). Способ построения отражен в английском названии графика. По осям графика откладываются значения квантилей: по оси OX – теоретических квантилей, квантилей эталонного нормального распределения; по оси OY – выборочных, эмпирических квантилей переменной, с которой мы работаем.

Для иллюстрации построения графика проведем предобработку данных. Давайте стандартизуем переменную роста – из каждого значения вычтем выборочное среднее и поделим на стандартное отклонение. В R это делается функцией scale(): \[height\_std=\frac{x_{i}-\bar{x}}{s_{x}}\]

С помощью этого преобразования мы добъемся того, что переменная будет принимать значения, похожие на значения нормального распределения.

df$height_std <- scale(df$height)
summary(df$height_std)
##        V1          
##  Min.   :-2.77698  
##  1st Qu.:-0.63748  
##  Median :-0.00141  
##  Mean   : 0.00000  
##  3rd Qu.: 0.63466  
##  Max.   : 2.42722

Так что же делает нормальная вероятностная бумага? Давайте найдем квантиль \(x_{0.1}\) для стандартного нормального распределения (если забыли, повторите, что такое квантиль):

qnorm(0.1, mean = 0, sd = 1)
## [1] -1.281552

А теперь давайте посмотрим, чему равен квантиль \(x_{0.1}\) для стандартизованной переменной роста:

quantile(df$height_std, 0.1)
##       10% 
## -1.273548

Первое значение, теоретический квантиль нормального распределения \(-1.281552\), будет отложено на оси OX. Второе значение, эмпирический квантиль для переменной \(-1.281552\), будет отложено на оси OY. Чем ближе эти значения друг к другу, тем ближе распределение переменной к нормальному! График q-q plot будет состоять из множества квантилей, которые в идеале должны выстроиться вдоль прямой \(y=x\).

Посмотрим на нормальную вероятностную бумагу для роста. Чтобы нарисовать график, нужна функция qqnorm(). Чтобы нарисовать прямую для удобства интерпретации, нужно реализовать команду qqline(), но только после того, как вы нарисовали сам график.

qqnorm(df$height_std)
qqline(df$height_std)

Мы действительно видим, что точки на графике расположены вдоль прямой. Наблюдается небольшая “лесенка”, когда при изменении значений по OX значение по OY не меняется. Это значит, что в выборке есть повторяющиеся значения: это вполне логично для переменной роста. Тем не менее, по графику видно, что значения повторяются по небольшому количеству раз.

Теперь давайте посмотрим на переменную, чье распределение не похоже на нормальное: количество баллов ЕГЭ по математике. Нарисуем гистограмму. Мы видим более массивный левый хвост, чем правый: более низких баллов больше, чем очень высоких. Кроме того, нет ярко выраженного “нормального пика”.

hist(df$math, freq = F, breaks = 25,
     col = "navajowhite", xlab = "Math score", main = "Histogram of math score")
curve(dnorm(x, mean(df$math, na.rm = T), sd(df$math, na.rm = T)), 
      add = T)

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

qqnorm(df$math)
qqline(df$math)

Статистическая оценка как случайная величина

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

Пердположим, сегодня 22 сентября, и москвичи и гости столицы бегут Московский Марафон. После того, как все участники финишируют, мы можем посчитать среднюю продолжительность забега для всех участников. Но спортивные результаты зависят от многих факторов помимо подготовки. Если мы возьмем всех тех же самых участников марафона и дадим им пробежать ту же самую дистанцию 23 сентября или 22 октября, то средняя продолжительность забега для всех участников уже будет другой. Это значит, что наша оценка средней продолжительности забега – случайная величина, и, если у нас много измерений, мы можем посчитать среднее, дисперсию, квантили…

Давайте посмотрим на примере.

На минутку представим, что участники нашего тренировочного опроса – единственные первокурсники-психологи на этой планете. Других нет. Легким движением руки наша база стала генеральной совокупностью.

Расчитаем среднее значение для роста респондентов в нашей базе. Будем считать, что это истинное среднее значение генеральной совокупности. Запомните его.

mean(df$height)
## [1] 169.0122

Теперь представим, что у нас нет данных по генеральной совокупности. Следовательно, мы не знаем и истинное среднее. В силу неизвестных причин мы не можем опросить всех первокурсников-психологов на этой планете, но можем опросить только 30 из них – и уже для них посчитать среднее значение роста. Это будет наша выборочная оценка истинного среднего.

sample1 <- sample(df$height, 30) # эта команда генерирует случайную выборку заданного объема 
mean(sample1)
## [1] 172.1667

Теперь сделаем следующее. Сгенерируем 1000 выборок размером 30 из нашей генеральной совокупности и для каждой из этих выборок посчитаем выборочную оценку роста. Так как выборки каждый раз разные, то и выборочные оценки будут разные: какие-то больше, какие-то меньше, чем истинное среднее.

# этот трюк выполняется каскадером, дома можно не повторять
# но если интересно, то эта штука называется loop 
means_sample <- c()
for (i in 1:1000){
  mean_sample <- mean(sample(df$height, 30))
  means_sample <- c(means_sample, mean_sample)
}

Посмотрим на описательные статистики для нашей 1000 выборочных средних. Минимум и максимум отличны от истинного среднего, но вот среднее выборочных средних стремится к истинному!

summary(means_sample)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   164.1   167.9   169.0   169.0   170.0   173.7

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

hist(means_sample, col = "plum", 
     freq= F, ylim = c(0,0.3), breaks = 15, xlab = "sample means of height", main = "Histogram of sample means")
curve(dnorm(x, mean(means_sample),
            sd(means_sample)), add = T)

Проверка статистических гипотез

Пример

Представьте, что вы живете в общежитии, и ваш сосед(ка) по комнате каждый раз, когда дело доходит до уборки, предлагает вам подбросить монетку. Если выпадает орел, убираете вы. Если выпадает решка, убирает он(а). Вы считаете, что это хорошее решение проблемы координации – случай слеп и непредвзят. Но к сотой уборке вы заметили, что 80 раз монетка выпадала орлом. Вы задумались: нужно ли в ваших бедствиях винить слепой случай или кривую монетку сосед(ки)?

Проверка статистических гипотез позволяет решить эту проблему.

Напомню, что в статистике мы живем в двух “параллельных” мирах. Первый мир – ненаблюдаемый – мир истинных распределений характеристик генеральной совокупности. Второй мир – наблюдаемый – мир выборочных оценок характеристик генеральной совокупности. Наши оценки что-то говорят нам об истинных распределениях, но насколько точно?

Но есть и хорошие новости. Пускай вы выдвинули некоторое предположение об истинном распределении. Пускай мы предположили, что монетка вашего(ей) соседа(ки) по общежитию правильная. Иными словами, орел и решка равновероятны. Это предположение называется нулевой гипотезой. \[H_{0}:~p=0.5\] Имеющиеся данные заставляют нас сомневаться в правдоподобности нулевой гипотезы. Кажется, орлы выпадают значительно чаще решек.

Мы можем, используя данные выборки, рассчитать некоторую характеристику этой выборки – статистику – которая будет количественно свидетельствовать о правдоподобности или неправдоподобности нулевой гипотезы. Мы знаем, что количество успехов в серии испытаний Бернулли описывается биномиальным распределением, а при достаточно больших n оно может быть приближено нормальным.

Пускай нашей статистикой будет стандартизованное значение количества выпавших орлов:

\[Z=\frac{80-n\times p}{\sqrt{npq}}=\frac{80-50}{5}=6\]

Мы получили значение случайной величины, которая имеет стандартное нормальное распределение. Обратите внимание, что значение нашей статистики – 6 – очень велико, оно находится далеко в правом хвосте стандартного нормального распределения. Вспомните, что полученное значение статистики, как и любая оценка, посчитанная по выборке, является случайной величиной. Если бы у нас была другая выборка, мы могли бы получить другое значение статистики. Но с какой вероятностью мы могли бы получить значение, большее данного, если мы считаем, что нулевая гипотеза верна? С какой вероятностью монетка бы обрекла вас на 95 из 100 дней уборки комнаты, если мы все еще думаем, что монетка правильная?

pnorm(80, mean = 50, sd = 5, lower.tail = F)
## [1] 9.865876e-10

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

Понятия

Мы уже знаем, что предположение об истинном распределении вероятностей называется нулевой гипотезой, \(H_{0}\). Нулевая гипотеза всегда простая и однозначно задает распределение. Мы всегда оцениваем правдоподобность нулевой гипотезы.

Помимо нулевой гипотезы, вам всегда требуется сформулировать альтернативную гипотезу, \(H_{a}\). Это конкурирующее предположение о распределении вероятностей. В случае с монеткой, нашей альтернативой было предположение \(H_{a} > 0.5\), что орлы выпадают чаще решек. Альтернативная гипотеза помогла нам рассчитать вероятность, которую мы использовали для заключения о монетке. Мы исключили ситуации, в которых монетка могла последовательно выпадать только решками. Такая альтернатива называется односторонней. Альтернативные гипотезы, которые имеют вид \(H_{a} \neq 0.5\), называются двухсторонними (как вы думаете, почему?).

Выборочная характеритика, позволяющая сделать вывод о правдоподобности гипотезы, называется статистикой критерия.

Статистический критерий – это инструмент, созданный для решения той или иной задачи: для проверки той или иной нулевой гипотезы. Поэтому вам всегда необходимо уметь сформулировать нулевую гипотезу. Критерии бывают параметрическими, использующимися, когда мы знаем распределение данных, и непараметрическими, использующимися, когда это распределение неизвестно.

Вероятность превысить наблюдаемое значение статистики, \(P(T>|t|)\) называется p-value. Также p-value определяется как вероятность отвергуть нулевую гипотезу, если она верна. Иными словами, с какой вероятностью мы ошибемся, если отвергем нулевую гипотезу?

Если p-value мало, то мы отвергаем нулевую гипотезу. Если p-value велико, то мы говорим, что мы не отвергаем нулевую гипотезу, или мы не имеем достаточных оснований, чтобы отвергуть нулевую гипотезу. Обратите внимание, что мы никогда не говорим, что мы принимаем нулевую гипотезу.

Проверка на нормальность

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

Критерий Колмогорова-Смирнова позволяет проверять гипотезу о равенстве выборочного распределения нормальному. Статистикой критерия выступает максимальное различие между эмпирической и теоретической функциями распределения.

Сгенерируем выброрку из стандартного нормального распредения с помощью функции rnorm(). В качестве аргумента указывается требуемый размер выборки.

set.seed(123) # эта команда позволяет воспроизводить результаты случайной генерации чисел
# если вы ее не прогоните, то сможете получить выборку, для которой критерий отвергнет нулевую гипотезу
sample2 <- rnorm(n = 200)

Проверим нашу выборку на нормальность. Мы видим, что p-value велик: у нас нет основания отвергать нулевую гипотезу о том, что выборка имеет нормальное распределение.

ks.test(sample2, "pnorm")
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sample2
## D = 0.053155, p-value = 0.6243
## alternative hypothesis: two-sided

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

ks.test(df$height, "pnorm", 
        mean = mean(df$height, na.rm = T), 
        sd = sd(df$height, na.rm = T))
## Warning in ks.test(df$height, "pnorm", mean = mean(df$height, na.rm = T), :
## ties should not be present for the Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  df$height
## D = 0.088671, p-value = 0.2882
## alternative hypothesis: two-sided

R выводит предупреждение, что присутствие повторяющихся значений может повлиять на рассчет p-value. Значение p-value велико, и мы не отвергаем нулевую гипотезу о нормальности.

Существует второй критерий проверки на нормальность: критерий Шапиро-Уилка.

shapiro.test(df$height)
## 
##  Shapiro-Wilk normality test
## 
## data:  df$height
## W = 0.98953, p-value = 0.475

P-value велико, и мы можем заключить, что рост респондентов распределен нормально.

Сравнение двух независимых групп

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

  1. t-критерий Стьюдента: параметрический критерий, используется, если есть нормальность распределения, и сравнивает средние в двух группах (используя значения переменных). \[H_{0}: a_{1}-a_{2}=0\]
  2. критерий Вилкоксона: непараметрический критерий, используется, если нет нормальности распределения, и сравнивает распределения в двух группах (используя ранги).

\[H_{0}: F(x)=G(x)\]

Следовательно, первым шагом в процедуре сравнения двух независимых групп будет проверка распределений по группам на нормальность.

Сравнение среднего роста девушек и юношей

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

library(ggplot2) # для удобства используем библиотеку ggplot2
ggplot(data = subset(df, gender != "NA"), # не учитываем респондентов, не указавших свой пол
       aes(y = height)) + geom_boxplot(fill = "plum") + facet_wrap(~gender)

Мы видим, что в выборке присутствуют нетипичные значения. Избавимся от них:

df_new <- subset(df, (height > 155 & gender == "male") | (height < 185 & gender == "female"))

Сначала проверим распределение роста на нормальность по группам. Визуально:

qqnorm(subset(df_new, gender == "female")$height)
qqline(subset(df_new, gender == "female")$height)

qqnorm(subset(df_new, gender == "male")$height)
qqline(subset(df_new, gender == "male")$height)

Распределение роста лежит примерно вдоль прямой.

Статистически:

shapiro.test(subset(df_new, gender == "female")$height)
## 
##  Shapiro-Wilk normality test
## 
## data:  subset(df_new, gender == "female")$height
## W = 0.98742, p-value = 0.4965
shapiro.test(subset(df_new, gender == "male")$height)
## 
##  Shapiro-Wilk normality test
## 
## data:  subset(df_new, gender == "male")$height
## W = 0.93884, p-value = 0.1873

Теперь нулевая гипотеза о нормальности не отвергается как для девушек, так и для юношей. Следовательно, мы будем использовать t-критерий Стьюдента.

Критерий Стьюдента

Есть две разновидности критерия Стьюдента:

  1. Дисперсии двух групп равны.
  2. Дисперсии двух групп неравны.

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

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

var.test(height ~ gender, # после знака ~ указывается группирующая переменная
         data = df_new, # база данных
         alternative = "two.sided" # вид альтернативы: выбираем двухстороннюю
         )
## 
##  F test to compare two variances
## 
## data:  height by gender
## F = 1.1245, num df = 95, denom df = 21, p-value = 0.7909
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5273772 2.0668978
## sample estimates:
## ratio of variances 
##            1.12453

Нулевая гипотеза не отвергается. Теперь мы можем воспользоваться критерием Стьюдента для равных дисперсий. Если бы на предыдущем шаге мы отвергли гипотезу о равных дисперсиях, мы бы указали аргумент var.equal = FALSE. Альтернативная гипотеза по дефолту двухсторонняя.

t.test(height ~ gender, df_new, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  height by gender
## t = -7.245, df = 116, p-value = 5.147e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -14.754199  -8.419097
## sample estimates:
## mean in group female   mean in group male 
##             167.1406             178.7273

p-value очень мало. Мы отвергаем нулевую гипотезу о равенстве средних между группами. Иными словами, молодые люди на первом курсе ОП “Психология”, кажется, в среднем выше девушек.

Критерий Вилкоксона

Если бы рост девушек и юношей не был бы распределен нормально, мы бы воспользовались критерием Вилкоксона:

wilcox.test(height ~ gender, df)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  height by gender
## W = 354.5, p-value = 3.83e-07
## alternative hypothesis: true location shift is not equal to 0