Математика и статистика, часть 2: семинар 15

Расчет p-value с помощью вероятностных функций

Напомню, что значение p-value, которое мы используем для принятия решения о нулевой гипотезе, представляет собой вероятность получить данное или еще более «нетипичное» значение статистистики критерия. Эта вероятность будет определяться по-разному в зависимости от того, работаете Вы с правосторонней альтернативной гипотезой: \[P(S\geq s_{\text{набл}})~,\] левосторонней альтернативной гипотезой: \[P(S\leq s_{\text{набл}})~,\] или двухсторонней альтернативной гипотезой: \[P(|S|\geq s_{\text{набл}})~,\] Таким образом, чтобы рассчитать p-value , нам нужно знать, какому распределению подчиняется статистика критерия, а также вид альтернативной гипотезы.

Давайте посмотрим, как осуществляется рассчет соответствующих вероятностей для известных нам распределений с помощью вероятностных функций R. Напомню, что функции вида pdist(), где вместо dist указывается необходимое распределение, возвращают значения функции распределения, или \(P(Z\leq z_i)\).

Например, для стандартноо нормального распределения (аргументы с параметрами можно не прописывать отдельно, так как по дефолту функция принимает именно их):

pnorm(1.96, mean = 0, sd = 1)

## [1] 0.9750021

Если же мы хотим рассчитать вероятность вида \(P(Z>z_i)\), то мы можем воспользоваться аргументом lower.tail, который позволяет «выбрать» нужный хвост распределения. Этот аргумент принимает логические значения: TRUE или FALSE (важно прописывать их заглавными буквами; также можно указывать только первую зглавную букву).

1-pnorm(1.96)

## [1] 0.0249979

pnorm(1.96, lower.tail = F)

## [1] 0.0249979

Если мы сложим две рассчитанные нами вероятности, то ожидаемо получим в сумме 1:

pnorm(1.96)+pnorm(1.96, lower.tail = F)

## [1] 1

Чтобы подсчитать \(P(|Z|\geq z_i)\), умножим вероятность верхнего хвоста на 2:

2*pnorm(1.96, lower.tail = F)

## [1] 0.04999579

Используя R, мы можем удобно рассчитывать аналогичные вероятности для распределения Стьюдента с помщью функции pt(). В качестве аргумента нам обязательно нужно указывать требуемое число степеней свободы. Для примера рассчитаем p-value для правосторонней альтернативы из задачки, разбиравшейся на семинаре, где \(t_{\text{набл}}=1.11\), а \(n=11\):

pt(1.11, df = 11-1, lower.tail = FALSE)

## [1] 0.1464885

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

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

Удобно посмотреть, как работает логика проверки гипотез, если мы будем знать истинное среднее. Давайте сгенерируем случайную выборку из стандартного нормального распределения. Так мы точно знаем, что \(\mu=0\).(Функция set.seed(7) здесь используется для воспроизводмости результатов генерации случайной выборки.)

set.seed(14)
s<-rnorm(100)

Рассчитаем получившееся выборочное среднее. Ожидаемо, оно не будет в точности равным 0:

mean(s)

## [1] 0.04314383

Можем ли мы считать, получив такую выборку, что истинное среднее равно 0? \[H_0: \mu = 0\] Давайте выберем правостороннюю альтернативу, так как выборочное среднее оказалось больше нашего гипотетического значения: \[H_a: \mu>0\] Давайте проверим сформулированную гипотезу с помощью критерия Стьюдента. В R это делается с помощью функции t.test(x = , mu = , alternative = ). В качестве аргументов требуется указать выборку x, гипотетическое значение mu и вид альтернативы alternative ("greater", "less" или "two.sided", все эти аргументы нужно указывать в кавычках):

t.test(s, mu = 0, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  s
## t = 0.47669, df = 99, p-value = 0.3173
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.1071319        Inf
## sample estimates:
##  mean of x 
## 0.04314383

В выдаче приводятся значения рассчитанной статистики критерия \(t = 0.47669\), количества степеней свободы \(df = n-1 = 99\) и p-value \(= 0.3173\). Мы получили достаточно большое значение p-value, которое подразумевает, что несмотря на то, что выборочное среднее отклонилось от нуля, наши данные не позволяют отвергнуть нулевую гипотезу. Можно предположить, что истинное среднее действительно равно 0 (а мы знаем это наверняка).

Обратите внимание, что функция t.test() также строит доверительный интервал, но делает это в соответствии с выбранной альтернативой. В нашем примере мы получили правосторонний доверительный интервал. Чтобы регулировать доверительную вероятность, надо воспользоваться аргументом conf.level внутри функции t.test().

Чип и Дейл спешат на помощь

Теперь попробуем проанализировать наши первые реальные данные. Мы провели опрос, где просили Вас поразмышлять о персонажах мультфильма «Чип и Дейл». Давайте посмотрим, что же у нас получилось.

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

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

Чтобы открыть файл с компьютера, можно сделать так:

df <- read.csv(file.choose())

В памяти R (environment) появился новый объект, который называется датафрейм (по простому — табличка с данными). Давайте на нее взглянем:

View(df)

Мы видим, что в нашей табличке 74 строки (74 респондента приняли участие в опросе) и 10 колонок (они обычно называются переменными).

Сейчас мы поработаем с Вашими ответами на вопрос, который звучал следующим образом: «Оцените свою ворчливость в процентах от Чипа (100% – такой же ворчливый, как Чип)». Давайте выведем описательные статистики для этой переменной. Чтобы выбрать колонку внутри таблицы, Вам нужно сначала указать название таблицы, затем прописать оператор $, как бы забираясь внутрь таблички, и далее указать название нужной колонки.

summary(df$chip)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   50.00   63.50   79.34   90.00 1000.00

Команда summary() приводит минимум, нижний квартиль, медиану, среднее, верхний квартиль и максимум нашей выборки. Если мы взглянем на максимум, то заметим, что он намного превышает «ожидаемый» максимум для этой переменной и очень сильно отклоняется от основного массива значений. Скорее всего, это выброс:

boxplot(df$chip, col = "plum")

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

df2 <- subset(df, chip < 1000)

Посмотрим на ящик с усами:

boxplot(df2$chip, col = "yellow1")

И на описательные статистики:

summary(df2$chip)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   50.00   62.00   66.72   90.00  101.00

Заметим, что выборочное среднее изменилось довольно сильно (было 79.34, а стало 66.72). Однако медиана почти не поменялась (было 63.5, стало 62). Действительно, выборочное среднее оказалось значительно менее устойчиво к нетипичным значениям, чем медиана (подумайте, почему).

Давайте попробуем проверить гипотезу о том, что студенты 1 курса ОП Политология в среднем ворчливы на 60% по сравнению с Чипом против правосторонней альтернативы (предполагая, что ворчливость в Чипах распределена нормально) на уровне значимости 5%:

t.test(x = df2$chip, mu = 60, alternative = "greater", conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  df2$chip
## t = 2.1605, df = 72, p-value = 0.01703
## alternative hypothesis: true mean is greater than 60
## 95 percent confidence interval:
##  61.53825      Inf
## sample estimates:
## mean of x 
##  66.72452

Получившийся p-value (или \(P(T\geq t_{\text{набл}}\)) равен 0.01703. Эта вероятность меньше, чем выбранный уровень значимости, поэтому мы отвергаем нулевую гипотезу. Судя по всему, ворчливость студентов выше, чем 60% от Чипа!

Попробуем сделать то же самое для двухсторонней альтернативы:

t.test(x = df2$chip, mu = 60, alternative = "two.sided", conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  df2$chip
## t = 2.1605, df = 72, p-value = 0.03406
## alternative hypothesis: true mean is not equal to 60
## 95 percent confidence interval:
##  60.51995 72.92910
## sample estimates:
## mean of x 
##  66.72452

Ожидаемо, p-value стал в два раза больше, так как сейчас он представлет собой вероятность \(P(|T|\geq t_{\text{набл}}\), и равен 0.03406. Но это значение тоже меньше уровня значимости, следовательно, мы опять отвергаем нулевую гипотезу.

Напомню, что мы считаем, что при \(n>30\) распределение Стьюдента довольно близко к стандартному нормальному распределению. Давайте сравним p-value из предыдущей выдачи, который был рассчитан с использованием t-распределения, и p-value, который мы могли бы получить, используя z-распределение. Наше значение статистики равно \(t_{\text{набл}}=2.1605\).

p_t <- 0.03406 # беру зачение p-value из предыдущей выдач
p_z <- round(2*pnorm(2.1605, lower.tail = F), 5) # округляю до 5 знака после запятой
p_z

## [1] 0.03073

p_t-p_z # сравним

## [1] 0.00333

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