Данный ноутбук является конспектом по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2019). Автор ноутбука вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
Ура! Мы добрались до самого главного. До проверки гипотез. Аналитики в разных компаниях проверяют их на ежедневной основе. Давайте попробуем разобраться зачем. Для начала немного кода для красивых картинок в ggplot2. Как обычно, вы не обязаны понимать этот код. Он тут нужен для красоты.
Как проверять гипотезы
Точные и асимптотические критерии
Гипотезы о долях
Гипотезы о средних
Гипотезы о дисперсиях
library("ggplot2") # Пакет для красивых графиков
library("grid") # Пакет для субплотов
# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
options(repr.plot.width=4, repr.plot.height=3)
# В этом блоке написаны несколько функций, которые рисуют красивые картинки для проверки гипотез.
# Можно не обращать на них внимание :)
# Область, которую надо будет закрасить на графике
limitRange <- function(fun, min, max) {
function(x) {
y <- fun(x)
y[x < min | x > max] <- NA
return(y)
}
}
z_stat_picture <- function(z_stat, alpha =0.05, alternative = 'two-sided'){
# Опции для размеров графика
options(repr.plot.width=6, repr.plot.height=3)
# Какие области надо закрасить (зависит от алтиернативы)
if(alternative == 'two-sided'){
z_crit <- qnorm(1-(alpha/2))
dlimit_left <- limitRange(dnorm, -Inf, -z_crit)
dlimit_right <- limitRange(dnorm, z_crit, Inf)
}
if(alternative == 'less'){
z_crit <- qnorm(alpha)
dlimit_left <- limitRange(dnorm, -Inf, z_crit)
dlimit_right <- limitRange(dnorm, Inf, Inf)
}
if(alternative == 'greater'){
z_crit <- qnorm(1 - alpha)
dlimit_left <- limitRange(dnorm, -Inf, -Inf)
dlimit_right <- limitRange(dnorm, z_crit, Inf)
}
# Основная картинка
p <- ggplot(data.frame(x=c(-4,4)), aes(x = x))+
stat_function(fun=dnorm) + # вся функция
stat_function(fun=dlimit_left, geom="area", fill="blue", alpha=0.2) +
stat_function(fun=dlimit_right, geom="area", fill="blue", alpha=0.2) +
geom_point(aes(x=z_stat,y=0), color="red", size=2) +
geom_vline(xintercept = z_stat, size=0.3, linetype="dashed", color="red") +
annotate("text", label=round(z_stat,2), x= z_stat, y=0.2, parse=T, size=4, color="darkred")
# Рисуем критические точки
if(alternative == 'two-sided'){
p <- p + annotate("text", label=round(z_crit,2), x=z_crit+0.1, y=-0.02, parse=T, size=4) +
annotate("text", label=round(-z_crit,2), x=-z_crit-0.1, y=-0.02, parse=T, size=4)
}
if(alternative == 'less'){
p <- p + annotate("text", label=round(z_crit,2), x=z_crit-0.1, y=-0.02, parse=T, size=4)
}
if(alternative == 'greater'){
p <- p + annotate("text", label=round(z_crit,2), x=z_crit+0.1, y=-0.02, parse=T, size=4)
}
# Немного заключительных настроек, связанных с темой и вывод на экран
p <- p + theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
return(p)
}
# Точно такой же рисунок для t-статистики
t_stat_picture <- function(t_stat, n, alpha =0.05, alternative = 'two-sided'){
# Опции для размеров графика
options(repr.plot.width=6, repr.plot.height=3)
# Какие области надо закрасить (зависит от алтиернативы)
if(alternative == 'two-sided'){
t_crit <- qt(1-(alpha/2),df=n)
dlimit_left <- limitRange(function(x) dt(x, df=n), -Inf, -t_crit)
dlimit_right <- limitRange(function(x) dt(x, df=n), t_crit, Inf)
}
if(alternative == 'less'){
t_crit <- qt(alpha,df=n)
dlimit_left <- limitRange(function(x) dt(x, df=n), -Inf, t_crit)
dlimit_right <- limitRange(function(x) dt(x, df=n), Inf, Inf)
}
if(alternative == 'greater'){
t_crit <- qt(1 - alpha,df=n)
dlimit_left <- limitRange(function(x) dt(x, df=n), -Inf, -Inf)
dlimit_right <- limitRange(function(x) dt(x, df=n), t_crit, Inf)
}
# Основная картинка
p <- ggplot(data.frame(x=c(-4, 4)), aes(x = x))+
stat_function(fun=dt, args = list(df=n)) + # вся функция
stat_function(fun=dlimit_left, geom="area", fill="blue", alpha=0.2) +
stat_function(fun=dlimit_right, geom="area", fill="blue", alpha=0.2) +
geom_point(aes(x=t_stat,y=0), color="red", size=2) +
geom_vline(xintercept = t_stat, size=0.3, linetype="dashed", color="red") +
annotate("text", label=round(t_stat,2), x= t_stat, y=0.2, parse=T, size=4, color="darkred")
# Рисуем критические точки
if(alternative == 'two-sided'){
p <- p + annotate("text", label=round(t_crit,2), x=t_crit+0.1, y=-0.02, parse=T, size=4) +
annotate("text", label=round(-t_crit,2), x=-t_crit-0.1, y=-0.02, parse=T, size=4)
}
if(alternative == 'less'){
p <- p + annotate("text", label=round(t_crit,2), x=t_crit-0.1, y=-0.02, parse=T, size=4)
}
if(alternative == 'greater'){
p <- p + annotate("text", label=round(t_crit,2), x=t_crit+0.1, y=-0.02, parse=T, size=4)
}
# Немного заключительных настроек, связанных с темой и вывод на экран
p <- p + theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
return(p)
}
# Точно такой же рисунок для хи-квадрат
chi_stat_picture <- function(chi_stat, n, alpha =0.05, alternative = 'two-sided'){
# Опции для размеров графика
options(repr.plot.width=6, repr.plot.height=3)
# Какие области надо закрасить (зависит от алтиернативы)
if(alternative == 'two-sided'){
chi_crit1 <- qchisq(1-(alpha/2), df=n)
chi_crit2 <- qchisq(alpha/2, df=n)
dlimit_left <- limitRange(function(x) dchisq(x, df=n), 0, chi_crit2)
dlimit_right <- limitRange(function(x) dchisq(x, df=n), chi_crit1, Inf)
}
if(alternative == 'less'){
chi_crit2 <- qchisq(alpha,df=n)
dlimit_left <- limitRange(function(x) dchisq(x, df=n), 0, chi_crit2)
dlimit_right <- limitRange(function(x) dchisq(x, df=n), Inf, Inf)
}
if(alternative == 'greater'){
chi_crit1 <- qchisq(1 - alpha,df=n)
dlimit_left <- limitRange(function(x) dchisq(x, df=n), 0, 0)
dlimit_right <- limitRange(function(x) dchisq(x, df=n), chi_crit1, Inf)
}
# Основная картинка
p <- ggplot(data.frame(x=c(0, qchisq(1-(alpha/2), df=n) + 10)), aes(x = x))+
stat_function(fun=dchisq, args = list(df=n)) + # вся функция
stat_function(fun=dlimit_left, geom="area", fill="blue", alpha=0.2) +
stat_function(fun=dlimit_right, geom="area", fill="blue", alpha=0.2) +
geom_point(aes(x=chi_stat,y=0), color="red", size=2) +
geom_vline(xintercept = chi_stat, size=0.3, linetype="dashed", color="red") +
annotate("text", label=round(chi_stat,2), x= chi_stat, y=0.2, parse=T, size=4, color="darkred")
# Рисуем критические точки
if(alternative == 'two-sided'){
p <- p + annotate("text", label=round(chi_crit1,2), x=chi_crit1+0.1, y=-0.02, parse=T, size=4) +
annotate("text", label=round(chi_crit2,2), x=chi_crit2-0.1, y=-0.02, parse=T, size=4)
}
if(alternative == 'less'){
p <- p + annotate("text", label=round(chi_crit2,2), x=chi_crit2-0.1, y=-0.02, parse=T, size=4)
}
if(alternative == 'greater'){
p <- p + annotate("text", label=round(chi_crit1,2), x=chi_crit1+0.1, y=-0.02, parse=T, size=4)
}
# Немного заключительных настроек, связанных с темой и вывод на экран
p <- p + theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
return(p)
}
На самом деле мы можем сказать, что статистическая гипотеза это абсолютно любое утверждение, которое возникает в нашей голове, если его можно проверить с помощью данных.
В Австралии женщин дискриминируют на рынке труда!
В Питере люди любят кофе больше, чем в москве!
Посетители постоянно покупают в паре пиво и подгузники. Мы могли бы использовать это, чтобы увеличить прибыль.
Осьминог Пауль - настоящий прорицатель! Его предсказания всегда сбываются!
Доходности по акциям Яндекса имеют нормальное распределение!
Каждое из этих утверждений может быть проверено на данных. Сразу же стоит сказать важную вещь. Гипотезы нельзя принмать! Если кто-то при вас сказал, что он принял гипотезу, нужно срочно ударить этого человека по губам. Предположим, что мы собрали данные, связанные с предсказаниями Пауля. Предположим, что данные не противоречат тому, что Пауль действительно видит будущее. Означает ли это, что Пауль - провидец?
Вовсе нет, это означает, что текущие данные не противоречат этому. Возможно, если мы ещё понаблюдаем за предсказаниями Пауля, гипотеза о том, что он провидец, отвергнется. Если данные не противоречат нашему утверждению, это не является доказательством его верности. В матемтике да это всегда да, а нет это хотя бы один раз нет. Возможно, наше нет ещё не наступило. Именно поэтому гипотезы нельзя принимать, но можно не отвергать. Ещё раз, ещё раз если кто-то говорит, что принял гипотезу, нужно уведомить его, что он её вовсе не принял, а не отверг. Не стесняйтесь делать это.
Можно собирать данные и смотреть не противоречат ли им гипотезы
Это не очень достоверно, надо быть готовым к ошибкам
Гипотезы нельзя принмать, их не отвергают
Давайте посмотрим на конкретный пример. Джеймс Бонд утверждает, что он предпочитает пить мартини взболтанным, но не смешанным. Чтобы проверить это на практике, можно предложить ему пройти так называемый blind test, или слепое тестирование. Можно было бы завязать ему глаза, несколько раз предложить на выбор взболтанный и смешанный мартини, а после этого спросить какой напиток он предпочитает В данном случае если бы Джеймс Бонд выбрал взболтанный напиток, это считалось бы успехом, потому что его выбор соответствовал бы его утверждению. В противном случае считалось бы, что произошла неудача.
Завязываем агенту глаза, предлагаем взболтанный и смешанный мартини
Если больше понравился взболтанный $X_i = 1$
Сформулируем нашу первую гипотезу. Обычно, гипотезу, выдвинутую для проверки называют нулевой или основной.
$$H_0: \quad \text{ Джеймс Бонд не различает два вида напитков и выбирает наугад.}$$На языке математики это утверждение можно сформулировать с посощью Бернулиевской случайной величины. Как мы помним, у неё есть единтвенный параметр $p$, вероятность успеха. Если Джеймс Бонд не различает два напитка и берёт их наугад, то $p = 0.5$. Если нет, то он должен выбирать правильный напиток чаще неправильного и $p > 0.5$. Нулевая гипотеза ровно так и будет выглядеть:
$$ H_0: \quad p = \frac{1}{2} $$Эту гипотезу мы должны проверять против какой-то альтернативы. Альтернатива может быть разной. С одной стороны, мы можем рассмотреть двустороннюю альтернативу, состояшую в том, что агент Бонд различает два напитка между собой
$$ H_1: \quad p \ne \frac{1}{2}. $$С другой, мы могли бы рассмотреть одну из односторонних альтернатив. Например, что Бонд предпочитает именно вболтанный Мартини
$$ H_1: \quad p > \frac{1}{2}. $$Выбока: $X_1, \ldots, X_n$. Если Бонд не различает два напитка, он берёт их наугад. Сформулируем это в терминах $X \sim Bern(p)$.
$$ \begin{aligned} H_0: & \quad p = 0.5 \\ H_1: & \quad p \ne 0.5 \end{aligned} $$Эксперимент с завязанными глазами нужно провести $n$ раз. В итоге у нас в руках окажется выборка, по которой мы сможем найти $\hat p$. Предположим, что на основе 10 испытаний мы получили, что $\hat p = \frac{6}{10}$. Означает ли это, что Бонд отвечает за свои слова? Возникает соблазн сказать, что да, так как число получилось больше $\frac{1}{2}$, и это ошибочный соблазн. Давайте вспомним немаловажный факт: наша оценка $\hat p$ это функция от выборки. При других наблюдениях мы могли бы получить другую оценку. Хочется понять насколько сильно эта оценка отличается от нашей.
Пусть $10$ экспериментов показали, что $\hat{p}=0.6$
Джеймс Бонд различает?
Да! Число 0.6 больше, чем 0.5!
Секундочку... Мы же уже умеем отвечать на этот вопрос. За точность оценки отвечает доверительный интервал. Мы помним, что
Исходя отсюда мы можем построить доверительный интервал и всё станет понятнее. Для нашей ситуации получим:
p_hat = 0.6
n = 10
sd_p = sqrt(p_hat*(1-p_hat)/n)
p_left = p_hat - qnorm(1-0.05/2)*sd_p
p_right = p_hat + qnorm(1-0.05/2)*sd_p
cat('Параметр p с вероятностью 95% лежит между',p_left, 'и', p_right, '\n')
cat('Ширина интервала:', p_right - p_left)
Параметр p с вероятностью 95% лежит между 0.2963637 и 0.9036363 Ширина интервала: 0.6072726
Диапозон вышел довольно широким. Более того, в него попадает $0.5$. Получается, что Джеймс вполне может врать.
Итак, для разных выборок, наша оценка способностей Бонда различать напитки, может оказаться разной. Если расстояние от оценки до предполагаемого значения, $\hat p - 0.5$, окажется довольно маленьким, судя по всему, Бонд не может распознавать напитки. Выбор между ними происходит абсолютно случайно.
Осталось только формализовать, когда расстояние от нашей оценки $\hat p$ до $0.5$ маленькое. Так как наша оценка - случайная величина, то и расстояние до $0.5$ будет слачейной величиной.
Мы знаем, что
$$ \hat p \overset{asy}{\sim} N \left(p, \frac{p(1-p)}{n} \right) $$Если нулевая гипотеза верна, тогда истиное значение $p$ будет равно $0.5$ и
$$ \hat p \overset{asy}{\sim} N \left(0.5, \frac{0.5(1-0.5)}{10} \right) $$Расстояние до $0.5$ будет распределено с нулевым математическим ожиданием
$$ \hat p - 0.5 \overset{asy}{\sim} N \left(0, \frac{0.5(1-0.5)}{10} \right) $$Поделим всё это для удобства на стандартное отклонение и получим статистику для проверки гипотезы
$$ \frac{\hat p - 0.5}{\sqrt{\frac{0.5(1-0.5)}{10}}} \overset{asy}{\sim} N (0,1). $$Мы выяснили как распределно "расстояние" до значения $0.5$ при верности нулевой гипотезы
Осталось выбрать порог, при котором мы будем считать, что $H_0$ верна
Если альтернатива двусторонняя, очень большие и очень маленькие значения статистики говорят против нашей гипотезы
На основе этой статистики мы можем судить о том, насколько вероятно было получить по выборке $\hat p$, если в реальности природа задумала $0.5$. При справедливости альтернативы более вероятны большие значения статистики. Дело осталось за малым - выбрать порог, с которого начинаются такие значения.
При проверке гипотезы мы можем ошибиться и отвергунть $H_0$, когда она была верна. Хотелось бы, чтобы количество таких ошибок было заранее зафиксировано. Будем перед проверкой гипотезы объявлять желаемое число таких ошибок - уровень значимости, $\alpha$.
Напомню, что если мы упаковываем $100$ парашютов на $5\%$ уровне значимости, значит $5$ человек умрут. Если мы зафиксировали для запаковки такой уровень значимости, значит мы готовы смириться с этими смертями.
Если альтернатива двусторонняя, то не хотелось бы получить очень маленькие значения статистики или, наоборот, очень большие её значения. Будем считать, что если наблюдаемое значение статистики, $z$ попало в $1-\alpha$ процентный доверительный интервал, то всё окей и данные не противоречат нулевой гипотизе.
Если мы отвергаем $H_0$, когда она верна - мы ошибаемся
Фиксируем число ошибок, $\alpha$ - уровень значимости (обычно берут 1% 5% или 10%)
Если мы упаковываем $100$ парашютов на уровне значимости $5\%$, то пять парашютистов умрут, и мы готовы с этим смириться
Если статистика попала в $1-\alpha$% доверительный интервал, то данные не противоречат гипотезе
Посмотрим на то как это выглядит на картинке. Мы зафиксировали $5\%$ уровень значимости. После посчитали наблюдаемое значение статистики. Оно оказалось равно $0.65$. Критические значения статистики равны $\pm 1.96$. Мы можем найти их исходя из нормального распределения, также как мы делали это для доверительных интервалов. За ними начинается критическая область. Если наблюдаемое значение попадает в эту область, то мы считаем расстояние между $\hat p$ и $0.5$ слишком большим.
p_hat = 0.6
n = 10
sd_p = sqrt(p_hat*(1-p_hat)/n)
z_stat = (p_hat - 0.5)/sd_p # считаем наблюдаемую статистику
z_stat_picture(z_stat, alpha = 0.05, alternative = 'two-sided')
Для нашей сиуации гипотеза о том, что Бонд не различает между собой напитки, не отвергается. Можно выбрать одностороннюю альтернативу
$$ H_1: \quad p > \frac{1}{2}. $$В таком случае против нулевой гипотезы будут свидетельствовать только большие значения статистики и мы все $\alpha\%$ критиической области должны перекинуть направо. Гипотеза о рандомности выбора Бонда снова не отвергается.
z_stat_picture(z_stat, alpha = 0.05, alternative = 'greater')
Голубая область - критическая область, при попадании статистики в неё - гипотеза отвергается, "расстояние" оказывается большим
В нашем случае наблюдаемое значение статистики попало в доверительную область
Гипотеза о том, что Бонд не различает между собой напитки, не отвергается для обеих альтернатив
Подведём итоги. Проверка гипотез заключается в понимании того, какая из двух альтернатив наиболее вероятно при имеющихся данных. Она осуществляется с помощью специально подобранных случайных величин, статистик. Статистики подбираются так, чтобы их распределение, при справделивости нулевой гипотезы, было известно. По сути, каждый критерий для проверки гипотезы это полноценная теорема со своими предпосылками, которые могут нарушаться. Алгоритм проверки гипотез выглядит следующим образом:
Проверяя гипотезы, мы можем совершать ошибки. Выделяют две ошибки: первого рода и второго рода. Мы можем либо ошибочно отвергнуть верную гипотезу, либо ошибочно не отвергнуть неверную гипотезу.
$\mbox{ }$ | $H_0$ верна | $H_0$ не верна |
---|---|---|
$H_0$ не отвергается | ok | Ошибка $2$ рода |
$H_0$ отвергается | Ошибка $1$ рода | ok |
Вероятность ошибки первого рода, $P(H_0 \text{ отвергнута} \mid H_0 \text{ верна})$, называют уровнем значимости и обозначают как $\alpha$.
Вероятность ошибки второго рода, $P(H_0 \text{ не отвергнута} \mid H_0 \text{ не верна})$, обозначают как $\beta$. Величину $1-\beta$ называют мощностью критерия. Критерий с наибольшей мощностью - классный.
Ошибки первого и второго рода неравнозначны, ошибка первого рода критичнее (фиксируем её в начале)
Ошибка второго рода минимизируется по остаточному принципу
Презумпция невиновности: подсудимый по умолчанию невиновен. Если нет доказательств обратному, нельзя утверждать, что он преступник, даже если он совершил преступление.
Ещё одно важное определение. P-значение или достигаемый уровень значимости облегчают проверку гипотез и выводятся почти во всех протоколах. В дальнейшем вы регулярно будете сталкиваться с этим понятием.
Сначала про то как это выглядит на пальцах. Выше мы с вами для двусторонней альтернативы построили красивую картинку. Давайте немного дополним её. Как обычно, код для рисования картинок можно не понимать.
dlimit_left <- limitRange(dnorm, -Inf, -z_stat)
dlimit_right <- limitRange(dnorm, z_stat, Inf)
z_stat_picture(z_stat, alpha = 0.05, alternative = 'two-sided') +
# Делаем красную заливку
stat_function(fun=dlimit_left, geom="area", fill="red", alpha=0.1) +
stat_function(fun=dlimit_right, geom="area", fill="red", alpha=0.1) +
# Добавляем симметричкую пунктирную линию
geom_point(aes(x=-z_stat,y=0), color="red", size=2) +
geom_vline(xintercept = -z_stat, size=0.3, linetype="dashed", color="red") +
annotate("text", label=round(-z_stat,2), x= -z_stat, y=0.2, parse=T, size=4, color="darkred")
Цифры $1.96$ и $-1.96$ это критические значения статистики. Голубая площадь за ними в сумме даёт уровень значимости $\alpha$. Красные точки это наблюдаемое значение статистики $0.65$ и симметричное ему значение $-0.65$. Красная плошадь это и есть P-значение.
Если красная площадь получилась больше синей, это означает, что наблюдаемое значение статистики попало между критическими и гипотеза не отвергается. Если красная площадь оказалась меньше синей, то это означает, что наблюдаемое значение ударило в один из хвостов и гипотеза отвергается. Посчитаем для данной ситуации p-значение.
z_stat
2*pnorm(-z_stat)
Получается, что красная площадь равна $0.52$. Это гораздо больше, чем $\alpha = 0.05$. Значит гипотеза не отвергается. Чем ниже получилось p-значение, тем сильнее данные свидетельствуют против нулевой гипотезы в пользу альтернативы.
p-значение упрощает проверку гипотез. В нашем случае оно больше уровня значимости, следовательно гипотеза не отвергается.
Вопрос: а какой уровень значимости надо выбрать, чтобы гипотеза впервые отверглась?
Ответ: равный p-значению.
Достигаемый уровень значимости или P-значение - это вероятность при справедливости нулевой гипотезы получить такое же значение статистики, как в эксперименте или ещё более экстремальное.
$$p = P(z \ge z_{\text{набл.}} \mid H_0)$$В нашем случае $p = 0.518$. То есть вероятность получить наши или ещё более экстримальные данные при условии справедливости нулевой гипотезы довольно высока. Это говорит в пользу нулевой гипотезы. Если бы вдруг получилось, что $p = 0.002$, то вероятность получить наши или ещё более экстримальные данные при верности нулевой гипотезы была бы очень маленькой и это говорило бы против нулевой гипотезы. Тем не менее, редкие события тоже происходят.
Обычно нас интересуют средние, доли и дисперсии
Будем проверять гипотезы о них
Каждый критерий для проверки гипотезы - небольшая теорема со своими предпосылками
Сегодня посмотрим параметрические критерии
В гипотезах для таких критериев высказывается предположение о значении параметра распределения, из которого предположительно поступила выборка
Выборка:
$$X_1, \ldots, X_n \sim iid \hspace{1mm} Ber(p)$$Нулевая и альтернативная гипотезы:
$$ \begin{aligned} &H_0: \hspace{2mm} p = p_0 \\ &H_1: \hspace{2mm} p < \ne > p_0 \end{aligned} $$Асимптотическое распределение оценки $\hat p = \bar x$ при верности нулевой гипотезы:
$$ \hat p \overset{asy}{\sim} N\left(p_0, \frac{p_0(1-p_0)}{n} \right) $$Статистика для проверки гипотезы:
$$ Z = \frac{\hat p - p_0}{\sqrt{\frac{p_0(1-p_0)}{n}}} \overset{asy}{\sim} N(0,1) $$В 70-х годах известный педиатр и автор книг по воспитанию детей Бенджамин Спок был арестован за участие в антивоенной демонстрации. Его дело должен был рассматривать суд присяжных.
Присяжные назначаются многоступенчатой процедурой. На очередном этапе остаётся 300 человек, из которых отбираются финальные 12.
В процессе Бенджамина Спока среди этих 300 только 90 были женщинами
Поскольку в те времена воспитанием детей занимались в основном женщины, Бенджамин Спок среди них был более популярен, поэтому адвокаты заподозрили, что обвинение специально пытается сделать финальный состав присяжных менее благосклонным к подсудимому.
Протест адвокатов. Обвинение специально пытается сделать итоговый состав благосклонным к подсудимому, убирая женщин.
Был ли отбор беспристрастным?
# Что делать?
Нулевая гипотеза состоит в том, что процедура отбора беспристрастна, женщины попадают в присяжные с вероятностью $0.5$. Подозрение состоит в том, что для женщин эта вероятность занижена. Из-за этого будем рассматривать в качестве альтернативной гипотезы - одностороннюю.
$$ \begin{aligned} &H_0: \hspace{2mm} p = 0.5 \hspace{2mm} \text{Отбор беспрестрастен} \\ &H_1: \hspace{2mm} p < 0.5 \hspace{2mm} \text{Женщин берут реже} \end{aligned} $$n = 300
p_hat = 90/300
z_stat = (p_hat - 0.5)/sqrt(0.5*(1-0.5)/n)
p_value = pnorm(z_stat)
cat('p-value:',p_value)
z_stat_picture(z_stat, alpha = 0.05, alternative = 'less')
p-value: 2.131096e-12
Видим, что наблюдаемое значение статистики оказалось глубоко в критической области. Гипотеза о беспрестрастности отбора отвергается. Значание $pvalue$ оказалось меньше $0.05$.
Можно попробовать использовать двустороннюю альтернативу. Это будет означать, что вероятность отбора женщин может быть как занижена, так и завышена.
n = 300
p_hat = 90/n
# Вручную
z_stat = (p_hat - 0.5)/sqrt(0.5*(1-0.5)/n)
z_stat_picture(z_stat, alpha = 0.05, alternative = 'two-sided')
p_value = 2*pnorm(z_stat)
cat('p-value:',p_value)
p-value: 4.262192e-12
Наши выводы это никак не поменяет. Сделаем то же самое во встроенном в R пакете. Это позволяеь сделать команда prop.test
.
test_res = prop.test(90, 300, p = 0.5, alternative = "less",
conf.level = 0.95, correct = FALSE)
test_res
1-sample proportions test without continuity correction data: 90 out of 300, null probability 0.5 X-squared = 48, df = 1, p-value = 2.131e-12 alternative hypothesis: true p is less than 0.5 95 percent confidence interval: 0.0000000 0.3451483 sample estimates: p 0.3
На выходе получаем табличку с основными результатами проверки гипотезы. В том числ оценку параметра и доверительный интервал для него. Через знак доллара можно вытянуть отдельные части протокола.
test_res$p.value
test_res$estimate
В случае двух выборок будем поступать ровно также, как при строительстве доверительных интервалов. Попробуйте проделать всё, что написано ниже самостоятельно руками.
Выборки:
$$ \begin{aligned} & X_1, \ldots, X_{n_1} \sim iid \hspace{1mm} Ber(p) \\ & Y_1, \ldots, Y_{n_2} \sim iid \hspace{1mm} Ber(p) \end{aligned} $$Нулевая и альтернативная гипотезы:
$$ \begin{aligned} &H_0: \hspace{2mm} p_1 = p_2 \\ &H_1: \hspace{2mm} p_1 < \ne > p_2 \end{aligned} $$Асимптотическое распределения оценок $\hat p_1 = \bar x$ и $\hat p_2 = \bar y$ при верности нулевой гипотезы:
$$ \hat p_1 \overset{asy}{\sim} N\left(p, \frac{p(1-p)}{n_1} \right) $$$$ \hat p_2 \overset{asy}{\sim} N\left(p,\frac{p(1-p)}{n_2} \right) $$Значит статистика $\hat p_1 - \hat p_2$ тоже будет иметь нормальное распределение. Найдём её параметры:
$$ \begin{aligned} & E(\hat p_1 - \hat p_2) = p - p = 0 & Var(\hat p_1 - \hat p_2) = p(1-p)\frac{n_1 + n_2}{n_1 \cdot n_2} \end{aligned} $$Статистика для проверки гипотезы:
$$ Z = \frac{\hat p_1 - \hat p_2}{\sqrt{P(1-P)\frac{n_1 + n_2}{n_1 n_2}}} \overset{asy}{\sim} N(0,1), $$где $P = \frac{m_1 + m_2}{n_1 + n_2} = \frac{\hat p_1 n_1 + \hat p_2 n_2}{n_1 + n_2}$
1600 гражданам Великобритании с правом голоса задают вопрос: одобряют ли они деятельность премьер-министра. 944 человека говорят, что одобряют. Через 6 месяцев опрос повторяется. На этот раз из 1600 опрошенных 880 говорят, что поддерживают премьер-министра. Нужно понять изменился ли рейтинг премьер-министра.
Чтобы понять, изменился ли рейтинг премьер-министра, нужно использовать статистический критерий. В качестве альтернативной гипотезы уместно взять двустороннюю, так как рейтинг мог как упасть, так и вырасти. Сначала проделаем всё вручную.
# Что делать?
В качестве альтернативы будет уместна двусторонняя гипотеза, так как рейтинг мог измениться в любую сторону.
$$ \begin{aligned} &H_0: \hspace{2mm} p_1 = p_2 \hspace{2mm} \text{Рейтинг не поменялся} \\ &H_1: \hspace{2mm} p_1 \ne p_2 \hspace{2mm} \text{Рейтинг поменялся} \end{aligned} $$n1 = 1600
m1 = 944
n2 = 1600
m2 = 880
p1_hat = m1/n1
p2_hat = m2/n2
P = (m1 + m2)/(n1 + n2)
z_stat = (p1_hat - p2_hat)/(sqrt(P*(1-P)*(n1+n2)/(n1*n2)))
z_stat_picture(z_stat, alpha = 0.05, alternative = 'two-sided')
p_value = 2*pnorm(-z_stat)
cat('p-value:',p_value)
p-value: 0.02229829
P-value меньше $0.05$. На уровне значимости $5\%$ гипотезе о том, что рейтинг премьера не изменился, отвергается. По-другому мы можем понять то же самое, увидем, что наблюдаемое значение статистики попадает в критическую область. Проверим гипотезу внутренними средствами R и получим тот же результат.
test_res = prop.test(x = c(m1, m2), n = c(n1, n2), correct = FALSE)
test_res
2-sample test for equality of proportions without continuity correction data: c(m1, m2) out of c(n1, n2) X-squared = 5.2224, df = 1, p-value = 0.0223 alternative hypothesis: two.sided 95 percent confidence interval: 0.005721636 0.074278364 sample estimates: prop 1 prop 2 0.59 0.55
В обоих опросах участвовали одни и те же люди, выборки связаны между собой
Z критерий предполагает, что выборки независимы
Нужен другой критерий
Мы можем перенести собранные по выборке данные в небольшую табличку
$\mbox{ } $ | + | - | $\Sigma$ |
---|---|---|---|
+ | 794 | 150 | 944 |
- | 86 | 570 | 656 |
$\Sigma$ | 880 | 720 | 1600 |
По столбцам отложим второй опрос, по строчкам первый. Из таблички можно понять, что $794$ человека поддержали премьера в обоих опросах, $150$ только в первом и т.д. Выше нам не требовалась такая конкретика, поэтому там мы не упоминал кто как голосовал.
$\mbox{ } $ | + | - | $\Sigma$ |
---|---|---|---|
+ | e | f | e+f |
- | g | h | g+h |
$\Sigma$ | e+g | f+h | n |
Если мы будем расчитывать итоговую статистику как
$$ Z = \frac{f - g}{\sqrt{f + g - \frac{(f-g)^2}{n}}} \overset{asy}{\sim} N(0,1), $$то мы внутри неё учтём только недиагональные элементы. То есть только те объекты, на которых значение в выборке различаются. Грубо говоря, это позволит нам подчистить нашу выборку от взаимосвязей.
n = 1600
f = 150
g = 86
z_stat = (f - g)/(sqrt(f + g - (f-g)^2/n))
z_stat_picture(z_stat, alpha = 0.05, alternative = 'two-sided')
p_value = 2*pnorm(-z_stat)
cat('p-value:',p_value)
p-value: 2.804011e-05
Гипотеза о том, что рейтинг премьера не изменился, после такой подчистки, начинает отвергаться более уверенно. Для проверки гипотезы о равенсве доле для связных выборок я не нашёл функции. Придётся реализовать её самостоятельно.
Проверять гипотезы о долях помогает Z-тест
Он имеет асимптотически нормальное распределение
Статистика для связных выборок и независимых считается немного по-разному
Обе имеют асимптотически нормальное распределение
Позволяет проверять гипотезы о математических ожиданиях нормального распределения. Предпосылка о наличии нормальности очень критична. Её нужно каждый раз проверять отдельно. На то, как это делается, мы посмотрим на следующей паре.
Выборка:
$$ X_1, \ldots, X_n \sim iid \hspace{1mm} N(\mu, \sigma^2) $$Нулевая и альтернативная гипотезы:
$$ \begin{aligned} &H_0: \hspace{2mm} \mu = \mu_0 \\ &H_1: \hspace{2mm} \mu < \ne > \mu_0 \end{aligned} $$Распределение оценки $\hat \mu = \bar x$ при верности нулевой гипотезы:
$$ \hat \mu \sim N\left(\mu_0, \frac{\sigma^2}{n} \right) $$Статистика для проверки гипотезы:
$$ Z = \frac{\hat \mu - \mu_0}{\frac{\sigma}{\sqrt{n}}} \sim N(0,1) $$Обратите внимание, что из-за нормальности распределения это распределение является точным, а не асимптотическим. Это очень важная деталь.
Обратите внимание на дисперсию. Про неё в нулевой гипотезе не делается никаких предположений. Она уже заранее должна быть для нас известна. Если это так, то с распределением всё хорошо и можно проверять гипотезу.
Что делать, если дисперсия неизвестна? Оценить её по выборке! Тогда у нас будет уже целых два распределения: для среднего и для выборочной дисперсии. Со средним всё просто:
С дисперсией чуток сложнее. По Тореме Фишера, если генеральная совокупность имеет нормальное распределение, то статистика
Не помните такого факта? Марш в лекции! Найдите его и внимательно изучите!
Владея этим фактом, мы можем сделать для статистики для проверки гипотезы небольшой финт, который позволит заменить неизвестное $\sigma$ на найденное по выборке $s$ и понять как будет распределена итоговая дробь:
Подытожим. Если дисперсия неизвестна и её надо оценить, то статистика для проверки гипотезы будет иметь распределение Стьюдента:
Статистика для проверки гипотезы:
$$ t = \frac{\hat \mu - \mu_0}{\frac{s}{\sqrt{n}}} \sim t(n-1) $$Обратите внимание, что это распределение точное. Также обратите внимание, что при $n \to \infty$ распределение Стьюдента сходится к нормальному. Значит наша статистика будет иметь асимптотически нормальное распределение и при больших $n$ можно смело использовать его для проверки гипотез. Погрешность при его использовании будет небольшой. Но ежели нужен точный результат, а не асимптотический, используйте распределение Стьюдента.
Средний вес детей при рождении 3300 г, у женщин, живущих за чертой бедности - 2800 г
25 женщин, живущих за чертой бедности, участвовали в экспериментальной программе ведения беременности
Средний вес их детей составил 3075 г, стандартное отклонение 500 г
Эффективна ли программа?
# Что делать?
В качестве альтернативы берём правостороннюю, так как вес детей должен был возрасти
$$ \begin{aligned} &H_0: \hspace{2mm} \mu = 2800 \hspace{2mm} \text{Программа не помогла}\\ &H_1: \hspace{2mm} \mu > 2800 \hspace{2mm} \text{Программа помогла} \end{aligned} $$Теоретическая дисперсия неизвестна. Её придётся оценивать. Из-за этого будем использовать $t-$статистику.
mu_0 = 2800
mu_hat = 3075
s = 500
n = 25
t_stat = sqrt(n)*(mu_hat - mu_0)/s
# обратите внимание, что p-value теперь мы считаем по
# распределению стьюдента! Картинку тоже строим для него
t_stat_picture(t_stat, n, alpha = 0.05, alternative = 'greater')
p_value = pt(-t_stat, df=n-1)
cat('p-value:',p_value)
p-value: 0.005573915
Гипотеза о том, что программа была неэффективной отвергается. Обратите внимание, что у $t$-распределения критические значения отличаются от нормального. Для данной ситуации это $1.71$, а не $1.64$. Напомню, что критическое значение $t$-статистики находится правее из-за более тяжёлых хвостов. Проверить то же самое внутренними средствами R помогает сделать команда t.test
. Чуть ниже мы увидим её применение.
Обратите внимание, при проверки гипотезы, мы неявно пользовались предположением, что вес детей распределён нормально. На практике это часто не так. В следующий раз мы научимся проверять гипотезу о нормальности выборки и ещё раз обсудим зачем она нужна.
Если выборка не распределена нормально, то наша оценка для среднего будет иметь асимптотически нормальное распределение:
$$ \hat \mu \overset{asy}{\sim} N\left(\mu_0, \frac{\sigma^2}{n} \right) $$В такой ситуации мы не можем легко нащупать точное распределение для статистики. Поэтому для проверки гипотезы мы будем использовать асимптотическое распределение статистики:
В табличке coal лежит статистика о авариях на угольных шахтах в Британии. Исследователь Андрей считает, что они происходят с интенсивностью $3$ штуки за год. Проверьте на $5\%$ уровне значимости гипотезу об этом.
# Загружаем табличку с данными
df = read.csv('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/coal.csv', sep=',')
head(df, 5) # шапка таблицы
X | year | count |
---|---|---|
1 | 1851 | 4 |
2 | 1852 | 5 |
3 | 1853 | 4 |
4 | 1854 | 1 |
5 | 1855 | 0 |
Гипотеза:
$$ \begin{aligned} &H_0: \hspace{2mm} \lambda = 3 \\ &H_1: \hspace{2mm} \lambda \ne 3 \end{aligned} $$Мы знаем, что $\hat \lambda = \bar x$. Воспользуемся асимптотическим критерием для средних. Как можно убедиться ниже, гипотеза отвергается.
n = nrow(df)
z_stat = (mean(df$count) - 3)/(var(df$count)/sqrt(n))
z_stat_picture(z_stat, alpha = 0.05, alternative = 'two-sided')
p_value = pt(z_stat, df=n-1)
cat('p-value:',p_value)
p-value: 8.811492e-07
Ещё раз, ещё раз! Критерий для среднего является точным только в случае нормального распределения! Обычно в матстате и метрике ищут нормальное распределение только для того, чтобы пользоваться точными распределениями. Больше ни для чего это не надо. Запомните это!
Двухвыборочные критерии Стьюдента позволяют сравнивать средние значение двух выборок из нормальных распределений. Вывести статистику и её разспеределение можно по аналогии с тем, как мы делали это в доверительных интервалах.
Выборки:
$$ \begin{aligned} & X_1, \ldots, X_{n_1} \sim iid \hspace{1mm} N(\mu_1, \sigma_1^2) \\ & Y_1, \ldots, Y_{n_2} \sim iid \hspace{1mm} N(\mu_2, \sigma_2^2) \end{aligned} $$Нулевая и альтернативная гипотезы:
$$ \begin{aligned} &H_0: \hspace{2mm} \mu_1 = \mu_2 \\ &H_1: \hspace{2mm} \mu_1 < \ne > \mu_2 \end{aligned} $$Асимптотическое распределения оценок $\hat \mu_1 = \bar x$ и $\hat \mu_2 = \bar y$ при верности нулевой гипотезы:
$$ \hat \mu_1 \sim N\left(\mu, \frac{\sigma_1}{n_1} \right) $$$$ \hat \mu_2 \sim N\left(\mu,\frac{\sigma_2}{n_2} \right) $$Значит статистика $\hat \mu_1 - \hat \mu_2$ тоже будет иметь нормальное распределение. Найдём её параметры:
$$ \begin{aligned} & E(\hat \mu_1 - \hat \mu_2) = \mu - \mu = 0 & Var(\hat \mu_1 - \hat \mu_2) = \frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2} \end{aligned} $$Статистика для проверки гипотезы:
$$ Z = \frac{\hat \mu_1 - \hat \mu_2}{\sqrt{ \frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}} } $$И вот тут начинаются проблемы.
Дисперсия считается по формуле:
$$ s^2 = \frac{(n_1-1) \cdot s_1^2 + (n_2-1) \cdot s_2^2}{n_1 + n_2 -2} $$Но есть две проблемы. Во-первых, число степеней свободы $\nu$ у этого распределения будет считатся по довольно сложной формуле:
$$ \nu = \frac{ \left(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2} \right)^2 }{ \frac{s_1^4}{n_1^2(n_1 -1)} + \frac{s_2^4}{n_2^2(n_2-1)}} $$Во-вторых, нулевое распределение будет не точным, а приближённым. Точного решения (то есть точного нулевого распределения для такой статистики) не существует. Это проблема называется проблемой Беренца-Фишера. Она говорит, что невозможно точно сравнить две средние величины для выборок, дисперсии которых неизвестны.
Апроксимация, которую мы посмотрели выше, будет достаточно точной в двух ситуациях: когда выборки одинаокового объёма. Либо если знак между $n_1$ и $n_2$ такой же , как между $\sigma_1$ и $\sigma_2$. То есть, если из выборки с большей дисперсией собрано больше наблюдений.
У вас в лекциях были первые две ситуации. Третью ситуацию вы не обсуждали, но я счёл, что вы должны о ней услышать.
Есть данные о продажной стоимости недвижимости в Сиэтле для 50 сделок в 2001 году и 50 в 2002. Изменились ли в среднем цены?
df = read.csv('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/seattle.csv', sep='\t')
head(df)
Price | Year |
---|---|
142.00 | 2002 |
232.00 | 2002 |
132.50 | 2002 |
200.00 | 2002 |
362.00 | 2002 |
244.95 | 2002 |
# Вытащим цены в отдельные векторы
price_2002 = df[df$Year == 2002,]$Price
price_2001 = df[df$Year == 2001,]$Price
# Расчитаем всё необходимое для проверки гипотезы
n_1 = length(price_2001)
n_2 = length(price_2002)
mu_1 = mean(price_2001)
mu_2 = mean(price_2002)
s_1 = var(price_2001)
s_2 = var(price_2002)
t_stat = (mu_1 - mu_2)/sqrt(s_1/n_1 + s_2/n_2)
Мы не можем сказать, что дисперсии известны или равны. Придётся использовать сложную статистику со сложной формулой. Либр придётся сначала проверить гипотезу о равенстве дисперсий, а уже затем проверять гипотезу о равенстве математических ожиданий.
a = (s_1/n_1 + s_2/n_2)^2
b = s_1^2/(n_1^2*(n_1-1)) + s_2^2/(n_2^2*(n_2-1))
a/b
t_stat_picture(t_stat, 72, alpha = 0.05, alternative = 'two-sided')
p_value = 2*pt(t_stat, df=72)
cat('p-value:',p_value)
p-value: 0.423054
Гипотза о равенстве цен не отвергается. Проверим эту же гипотезу внутренней командой R. Аргумент var.equal
отвечает за равенство дисперсий. Если нам откуда-то известно, что она равны, гипотеза будет проверяться с учётом этого. Можете на досуге попробовать.
t.test(price_2001, price_2002, var.equal = FALSE,
alternative = "two.sided", conf.level = 0.95)
Welch Two Sample t-test data: price_2001 and price_2002 t = -0.80573, df = 71.895, p-value = 0.4231 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -140.11588 59.45468 sample estimates: mean of x mean of y 288.9265 329.2571
Внутри протокола получаем точно такие же значения статистики, степеней свободы и p-значение, как при расчёте руками.
Измерения делаются на одних и тех же объектах
Это значит, что $\bar y - \bar x$ это то же самое, что $\overline{x-y}$
Благодаря тому, что мы можем смотреть прирост по каждому объекту, мы скатываемся к одновыборочному критерию
Если выборки связаны между собой, то ситуация очень сильно упрощается. Фактически это означает, что у нас был какой-то объект. Мы пронаблюдали по нему выборку $X_1, \ldots, X_n$. Потом мы что-то сделали с этим объектом и после пронаблюдаль выбоку $Y_1, \ldots, Y_n$. После мы решили проверить гипотезу о том, что наше воздействие имело силу.
В данном случае разность $\bar y - \bar x$ это то же самое, что $\overline{x-y}$. Таким образом, $t$-критерий для двух связных выборок эквивалентен одновыборочному $t$-критерий, примененному к выборке попарных разностей.
Из-за того, что мы наблюдаем одни и те же объекты, никаких проблем с дисперсией не возникает.
# Что делать?
Препараты принимали одни и те же дети. Выборка является связанной. Поэто будем использовать критерий стьюдента для связанных выборок.
$$ \begin{aligned} &H_0: \hspace{2mm} \mu_1 = \mu_2 \hspace{2mm} \text{Препарат неэффективен} \\ &H_1: \hspace{2mm} \mu_1 > \mu_2 \hspace{2mm} \text{Препарат хорош} \end{aligned} $$df = read.csv('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/ADHD.csv', sep=' ')
head(df)
D0 | D60 |
---|---|
57 | 62 |
27 | 49 |
32 | 30 |
31 | 34 |
34 | 38 |
38 | 36 |
d = df$D0 - df$D60
mu_d = mean(d)
s_d = sd(d)
n = nrow(df)
t_stat = mu_d/(s_d/sqrt(n))
t_stat_picture(t_stat, n-1, alpha = 0.05, alternative = 'two-sided')
p_value = 2*pt(t_stat, df=n-1)
cat('p-value:',p_value)
p-value: 0.003771488
Гипотеза о том, что лекарство было неэффективно, отвергается. Попробуем внутреннюю команду. Для этого нам понадобится ещё один аргумент, pared
, отвечающий за связанность выборок.
t.test(df$D0, df$D60, var.equal = FALSE, paired=TRUE,
alternative = "two.sided", conf.level = 0.95)
Paired t-test data: df$D0 and df$D60 t = -3.2224, df = 23, p-value = 0.003771 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -8.141431 -1.775236 sample estimates: mean of the differences -4.958333
Обратите внимание, что если не учесть, что выборки связанные, можно получить неправильный результат и сделать вывод, что лекарство не помогает детям.
t.test(df$D0, df$D60, var.equal = FALSE,
alternative = "two.sided", conf.level = 0.95)
Welch Two Sample t-test data: df$D0 and df$D60 t = -1.4522, df = 45.671, p-value = 0.1533 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -11.832593 1.915927 sample estimates: mean of x mean of y 39.75000 44.70833
Для проверки гипотезы о равенстве средних для произвольных распределений используют асимптотический тест:
$$ Z = \frac{\hat \mu_1 - \hat \mu_2}{\sqrt{ \frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}} } \overset{asy}{\sim} N(0,1) $$Если выборка распределена нормально, можно использовать t-статистику, имеющую распределение стьюдента
Если средних два и выборки независимы, возможны три разные ситуации
Если выборки связанные, то всё скатывается к одновыборочному критерию
Если выборка распределена не нормально, то используются те же статистики, но они имеют асимптотичеси нормальное распределение
Выборка: $X_1, \ldots, X_n \sim iid \hspace{1mm} N(\mu, \sigma^2)$
Нулевая и альтернативная гипотезы:
$$ \begin{aligned} &H_0: \hspace{2mm} \sigma^2 = \sigma^2_0 \\ &H_1: \hspace{2mm} \sigma^2 < \ne > \sigma^2_0 \end{aligned} $$Говоря другими словами, мы хотим понять насколько сильно дисперсия отличается от $\sigma_0^2$. Для того, чтобы измерить расстояние от полученной по выборке дисперсии до дисперсии, предположенной в нулевой гипотезе, нужна статистика.
Вспоминаем, о чём мы говорили в доверительных интервалах. Нужно придумать хорошего союзника. Для нормальных распределений хорошим оказывается теорема Фишера. Из-за неё обычно смотрят на отношение дисперсий:
$$\frac{s^2}{\sigma_0^2}.$$Если отношение оказалось близким к единице, то супер! Дисперсии очень похожи и гипотеза не отвергается. Если отношение слишком маленькое или слишком большое, то плохо. Несостыкуется. Как обычно, осталось придумать статистику, которая помогла бы нам понимать когда расстояние получилось слишком большим, а когда слишком маленьким.
Если математическое ожидание неизестно, то используем теорему Фишера и понимаем, что
$$\frac{(n-1) s^2}{\sigma^2} =\sum_{i=1}^n \left(\frac{X_i - \bar x}{\sigma}\right)^2 \sim \chi^2_{n-1}.$$Как вы помните из тетрадки про доверительные интервалы, одна степень свободы сгорает из-за того, что мы вынуждены оценивать сразу два параметра: среднее и дисперсию. Это накладывает на нашу выборку одно линейное ограничение и количество случайных величин, доступных для варьирования, уменьшается на одну. Обратите внимание, что из-за нормальности распределения эта статистика является точной, а не асимптотической. Это очень важная деталь.
Давайте ссимулируем выборку из нормального распределения и проверим обоими критериямигипотезу о дисперсии для него. Будем проверять гипотезу о том, что дисперсия равна $4$ против гипотезы о том, что нет.
y <- rnorm(10, mean = 10, sd = 2)
y[1:10]
Преположим, что нам откуда-то известно математическое ожидание.
n = length(y)
s2 = 1/n*sum((y-10)^2)
chi_stat = n*s2/4
chi_crit = qchisq(1-0.05/2,df=n)
chi_stat
chi_stat_picture(chi_stat, n)
Гипотеза о равенстве дисперсии $4$ не отвергается.
Теперь сделаем то же самое для случая, когда математическое ожидание неизвестно.
n = length(y)
s2 = var(y)
chi_stat = (n-1)*s2/4
chi_crit = qchisq(1-0.05/2,df=n-1)
chi_stat
chi_stat_picture(chi_stat, n-1)
Снова гипотеза не отвергается, но значения статистик при этом получаются другими.
Выборки:
$$ \begin{aligned} & X_1, \ldots, X_{n_1} \sim iid \hspace{1mm} N(\mu_1, \sigma_1^2) \\ & Y_1, \ldots, Y_{n_2} \sim iid \hspace{1mm} N(\mu_2, \sigma_2^2) \end{aligned} $$Нулевая и альтернативная гипотезы:
$$ \begin{aligned} &H_0: \hspace{2mm} \sigma^2_1 = \sigma^2_2 \\ &H_1: \hspace{2mm} \sigma^2_1 < \ne > \sigma^2_2 \end{aligned} $$По аналогии с предыдущим пунктом, будем смотреть на то, насколько велико расстояние между двумя дисперсиями:
$$ \frac{s_1^2}{s_2^2} $$Если эта величина окажется близка к единице, значит дисперсии похожи. Осталось найти статистику, которую можно было бы использовать, чтобы понять насколько выборочные дисперсии оказались близки. При верности нулевой гипотезы, обе выборки обладают одинаковой дисперсией, значит мы можем сделать следущий финт:
При верности нулевой гипотезы, получаем, что
$$ F \sim \frac{\chi^2_{n_1-1}}{n_1-1} : \frac{\chi^2_{n_2-1}}{n_2-1} = F(n_1-1,n_2-1). $$Отношение дисперсий, при верности нулевой гипотезы, имеет распределение Фишера с $n_1-1$ и $n_2-1$ степенями свободы. Именно его мы и будем использовать как статистику для проверки гипотезы. Ещё раз обращаю ваше внимание на то, что такое распределение будет только при нормальности выборок.
Вашему вниманию представлены результаты прыжков в длину Васи Сидорова на двух соревнованиях.
На первых среди болельщиц присутствовала Аня Иванова (его первая любовь).
На вторых Аня среди болельщиц не присутствовала.
Проверьте гипотезу, состоящую в том, что отсуствие Ани позволяет Васе лучше сконцентрировтаься и прыгать с большей точность.
# Что делать?
Пара гипотез:
$$ \begin{aligned} &H_0: \hspace{2mm} \sigma_1^2 = \sigma_2^2 \\ &H_1: \hspace{2mm} \sigma_1^2 < \sigma_2^2 \end{aligned} $$first = c(1.83, 1.64, 2.27, 1.78, 1.89, 2.33, 1.61, 2.31)
second = c(2.05, 1.07, 1.59, 1.96, 1.29, 1.52, 1.18, 1.47)
f_stat = var(first)/var(second)
f_crit = qf(0.05, df1 = length(first) - 1, df2 = length(second) - 1)
f_stat
f_crit
Видим, что гипотеза о том, что Вася прыгает на обоих соревнованиях с одинаковой точностью не отвергается. Аня не оказывает на него никакого влияния. То же самое можно вроверить внутри встроенной команды.
В R проверить гипотезы, связанные с равенством дисперсий нормального распределения, позволяет команда var.test
.
var.test(first, second, alternative = 'less', conf.level = 0.95)
F test to compare two variances data: first and second F = 0.74409, num df = 7, denom df = 7, p-value = 0.3532 alternative hypothesis: true ratio of variances is less than 1 95 percent confidence interval: 0.000000 2.817898 sample estimates: ratio of variances 0.7440892
Для нормальных выборок можно использовать статистики Фишера
Для двухвыборочного критерия подойдёт F-критерий
Есть и другие статитсики для проверки гипотезы о дисперсиях, но мы в них углубляться не будем
Интерес представляет не величина p-значения, а размер эффекта - степень отклонения данных от нулевой гипотезы
На основе размера эффекта делается вывод о том, являются ли результаты практически существенными
На самом деле эксперименты проводятся не для того, чтобы получить значение p-value. Как правило, исследователя интересует размер эффекта, то есть степень отклонения данных от нулевой гипотезы. Например, если эксперимент связан с проверкой способностей предсказателя будущего, то размер эффекта — это вероятность верного предсказания. Если проверяется эффективность лекарства, то размер эффекта — это вероятность выздоровления пациента, который это лекарство принимает, за вычетом эффекта плацебо. При запуске программы лояльности для пользователей интернет-магазина размер эффекта — это последующее увеличение среднего чека.
Размер эффекта — это величина, определенная на генеральной совокупности. Но, как правило, у исследователя есть только небольшая выборка из нее, а оценка размера эффекта по выборке — это случайная величина. Маленький достигаемый уровень значимости является показателем того, что такую оценку размера эффекта, какая получена по выборке, с маленькой вероятностью можно было получить случайно.
Достигаемый уровень значимости зависит не только от размера эффекта, но и от объема выборки, по которой оценивается эффект. Если выборка небольшая, скорее всего, нулевая гипотеза на ней не отвергается (если только она не слишком дикая). Однако с ростом объема выборки начинают проявляться все более тонкие отклонения данных от нулевой гипотезы. Велика вероятность, что на достаточно большой выборке значительная часть разумных нулевых гипотез будет отвергнута. Именно поэтому, даже если нулевая гипотеза отвергнута, это еще не значит, что полученный эффект имеет какую-то практическую значимость, её нужно оценивать отдельно. Чтобы лучше это понять, давайте рассмотрим несколько примеров.
за три года женщины, упражнявшиеся не меньше часа в день, набрали значимо меньше веса, чем женщины, упражнявшиеся меньше 20 минут в день $(p < 0.001)$
разница в набранном весе составила 150 г.
практическая значимость такого эффекта сомнительна
Первый пример связан с большим исследованием, в рамках которого на протяжении трех лет у большой выборки женщин измеряли вес, а также оценивали, насколько активно они занимаются спортом. По итогам исследования выяснилось, что женщины, которые в течение этого времени упражнялись не меньше часа в день, набрали значительно меньше веса, чем женщины, которые упражнялись менее $20$ минут в день. Статистическая значимость этого результата достаточно высока: $p < 0.001$. Проблема в размере эффекта: разница в набранном весе между двумя исследуемыми группами женщин составила всего $150$ граммов. $150$ граммов за $3$ года — это не очень много. Крайне сомнительно, что этот эффект имеет какую-то практическую значимость.
в 2002 году клинические испытания гормонального препарата Премарин, облегчающего симптомы менопаузы, были досрочно прерваны
было обнаружено, что его приём ведёт к значимому увеличению риска развития рака груди на $0.08\%$, риска инсульта на $0.08\%$ и инфаркта на $0.07\%$
формально эффект крайне мал, но с учётом численности населения он превращается в тысячи дополнительных смертей
Еще один пример связан с клиническими испытаниями гормонального препарата «Премарин» (в книге 2 глава, нет я не читал), который облегчает симптомы менопаузы. В $2002$ году эти испытания были прерваны досрочно, поскольку было обна- ружено, что прием препарата ведет к значимому увеличению риска развития рака груди (на $0.08\%$), инсульта (на $0.08\%$) и инфаркта (на $0.07\%$). Этот эффект статистически значим; при этом на первый взгляд кажется, что размеры эффектов ничтожны. Например, если кому-то сказать, что его любимые конфеты повышают риск возникновения инфаркта на $0.07\%$, вряд ли это заставит человека отказаться от этих конфет. Тем не менее, если пересчитать размеры эффектов на всю популяцию людей, которым этот препарат может быть потенциально приписан, результатом будут тысячи дополнительных смертей. Разработчики препарата не могут взять на себя эту ответственность, поэтому такой препарат немедленно запрещают и снимают с рынка. Этот пример показывает, что практическую значимость результата нельзя определить на глаз. В идеале она должна определяться человеком, который поставил задачу и понимает предметную область.
лекарство, замедляющее ослабление интеллекта больных Альцгеймером
при испытаниях выяснилось, что разница в IQ контрольной и тестовой групп составляет 13 пунктов
разница статистически нещначима
возможно, выборка мала для детекции эффекта и изучение лекарства стоит продолжить
Еще один пример — это испытание лекарства, которое замедляет ослабление интеллекта у людей, страдающих болезнью Альцгеймера. В этом исследовании очень сложно измерить размер эффекта. В течение эксперимента одна часть испытуемых должна принимать лекарство, а другая — плацебо. Только по прошествии нескольких лет можно будет сравнить эти две группы. Поэтому такое исследование длится долгое время и дорого стоит. Если при испытании оказывается, что разница между снижением IQ в контрольной группе, где люди принимали плацебо, и тестовой группе, где люди принимали препарат, составляет 13 пунктов, это различие очень большое, и на практике этот эффект крайне значим. При этом может оказаться, что ста- тистическая значимость не была достигнута, то есть $p > \alpha$, и формально нулевую гипотезу об отсутствии эффекта лекарства нельзя отвергнуть. Если предмет исследования очень важен, то, оказавшись в подобной ситуации, возможно, стоит продолжать исследования: набрать еще выборку, уменьшить дисперсию оценки размера эффекта и убедиться в том, что важное открытие не упущено.