Данный ноутбук является конспектом по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2019). Автор ноутбука вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
В прошлой тетрадке мы поговорили о том, как можно проверять гипотезы о долях, средних и дисперсиях с помощью параметрических тестов. В этой тетрадке мы попробуем обсудить ещё несколько важных моментов, сопряжённых с тестами. План будет следующим:
library("ggplot2") # Пакет для красивых графиков
library("grid") # Пакет для субплотов
library('rio')
# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
options(repr.plot.width=4, repr.plot.height=3)
Про ошибки 1 и 2 рода
Ещё разок про тест отношения правдоподобий
Бутстрэп
Множественная проверка гипотез
Параметрические критерии
Непараметрические критерии
Тест отношения правдоподобий
Бутстрэп
Мы уже говорили про тест отношения правдоподобий в прошлый раз, но вы все были в этот момент опустошены, поэтому давайте обсудим его ещё разок.
На них мы решили несколько задачек. Смотрите рукописные конспекты :)
Шаг первый: оценили модель методом максимального правдоподобия без ограничений и получили $L(\hat \theta_{ML})$
Шаг второй: оценили модель, наложив ограничение и получили $L(\hat \theta_{0})$
Оказывается, что
$$ LR = 2\cdot(\ln L(\hat \theta_{ML}) - \ln L(\theta_0)) \sim \chi^2_q, $$где $q$ — число ограничений
Пусть $X_1, \ldots, X_{30} \sim iid \hspace{1mm} Pois(\lambda)$ — количество серий ИП, которое смотрит Ульяна каждый день. Получилось, что $\sum X_i = 100$. Нужно проверить гипотезу о том, что
$$ H_0: \lambda = 2 \\ H_A: \lambda \ne 2 $$lnL <- function(l){
return(100*log(l) - 30*l)
}
2*(lnL(10/3) - lnL(2)) # наблюдаемое значение
qchisq(0.95, df=1) # критическое значение
df = import('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/flat.csv')
head(df)
n | price | totsp | livesp | kitsp | dist | metrdist | walk | brick | floor | code |
---|---|---|---|---|---|---|---|---|---|---|
1 | 81 | 58 | 40 | 6 | 12.5 | 7 | 1 | 1 | 1 | 3 |
2 | 75 | 44 | 28 | 6 | 13.5 | 7 | 1 | 0 | 1 | 6 |
3 | 128 | 70 | 42 | 6 | 14.5 | 3 | 1 | 1 | 1 | 3 |
4 | 95 | 61 | 37 | 6 | 13.5 | 7 | 1 | 0 | 1 | 1 |
5 | 330 | 104 | 60 | 11 | 10.5 | 7 | 0 | 1 | 1 | 3 |
6 | 137 | 76 | 50 | 9 | 11.0 | 7 | 1 | 1 | 1 | 8 |
x = log(df$price)
lnL <- function(param,data) {
n <- length(data)
result <- -n/2*log(2*pi)-n/2*log(param[2]^2)-sum((data-param[1])^2)/(2*param[2]^2)
return(result)
}
library('maxLik')
res <- maxLik(lnL, start=c(0.1,0.1),data=x)
summary(res)
-------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 16 iterations Return code 1: gradient close to zero Log-Likelihood: -539.0942 2 free parameters Estimates: Estimate Std. error t value Pr(> t) [1,] 4.791375 0.006978 686.66 <2e-16 *** [2,] 0.315159 0.004934 63.87 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 --------------------------------------------
# вбили функцию правдоподобия с ограничением, нашли её оптимум
lnL_res <- function(sigma,data) lnL(c(5, sigma),x)
res_restr <- maxLik(lnL_res, start=0.1,data=x)
summary(res_restr)
-------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 9 iterations Return code 2: successive function values within tolerance limit Log-Likelihood: -909.7547 1 free parameters Estimates: Estimate Std. error t value Pr(> t) [1,] 0.377955 0.005917 63.87 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 --------------------------------------------
2*(lnL(res$estimate, x) - lnL_res(res_restr$estimate, x))
qchisq(0.95, df=1)
А теперь давайте проверим гипотезу о том, что $\mu = 4.8$ и $\sigma^2 = 0.3$.
mu0 = 4.8
s0 = 0.3
2*(lnL(res$estimate, x) - lnL_res(c(mu0, s0), x))
qchisq(0.95, df=2)
В примерах выше мы сравнивали $\frac{L_{ML}}{L_R}$ с единицей
На основе этого мы делали вывод насколько сильно различаются два правдоподобия и решали что делать с гипотезой
Сравнивать отношение правдоподобий можно не с единицей, а с какой-то константой $c$
Лемма Неймана-Пирсона утверждает, что можно подобрать $c$ таким, чтобы у критерий была максимальная мощность (это было на лекции)
Подробнее смотри в Черновой страницы 86-67
Бутстрап - это метод оценки статистик сложных распределений. Часто возникает необходимость проверить гипотезу о какой-то очень неудобной статистике, распределение которой неизвестно или построить для неё доверительный интервал.
Чтобы сделать это, надо придумать как распределение этой статистики получить. Мы с вами уже знаем два способа сделать это.
Способ номер один называется параметрическим. Мы предполагаем, что генеральная совокупность имеет какое-то распределение. Мы можем даже проверить гипотезу об этом с помощью какого-нибудь критерия. После мы на основе этого распределения можем придумать адекватную случайную величину для проверки гипотезы. Именно так мы поступали раньше. Иногда такую случайную величину придумать сложно или вовсе невозможно. Например, поди придумай адекватный критерий для медианы или моды. К тому же не очень понятно из каких соображений выбирать семейство распределений для проверки гипотезы, так как про данные ничего неизвестно. Другим способом был тест отношения правдоподобий. Пришло время познакомиться с третьим методом.
Он говорит следующее: давайте извлечём из генеральной совокупности какое-то количество выборок, посчитаем по ним нашу статистику, для которой доверительный интервалы мы хотим получить, а потом оценим эмпирически её выборочное распределение и на его основе получим доверительные интервалы. Этот способ применим скорее в теории, чем на практике. Если мы можем генерировать бесконечное число выборок из генеральной совокупности, то для нас незатруднительно посчитать истинное значение интересующих нас параметров.
Эти мысли подводят нас к идее бутстрапа. В нашем распоряжении есть выборка. Давайте сделаем вид, что она и есть генеральная совокупность и будем извлекать из неё с повторением элементы. На основе получившихся подвыборок мы можем оценить всё, что нашей душе угодно. Поначалу такое кажется безумием, но это реально работает. И есть даже несколько теорем, которые доказывают почему это работает и как нужно правильно делать это в сложных ситуациях.
Опишем бутстрап чуть более формально. Пусть имеется выборка $X$ размера $N$. Равномерно возьмём из выборки $N$ объектов с возвращением. Это означает, что мы будем $N$ раз выбирать произвольный объект выборки (считаем, что каждый объект достаётся с одинаковой вероятность $\frac{1}{N}$), причём каждый раз мы выбираем из всех исходных $N$ объектов. Можно представить себе мешок, из которого достают шарики: выбранный на каком-то шаге шарик возвращается обратно в мешок, и следующий выбор опять делается равновероятно из того же числа шариков. Конечно же из-за возвращения среди них окажутся повторы.
Обозначим новую выборку через $X_1$. Повторяя процедуру $B$ раз, сгенерируем $B$ подвыборок $X_1, \ldots, X_B$. Теперь мы имеем достаточно большое число выборок и можем оценивать различные статистики исходного распределения.
По изначальной выборке мы могли посчитать всего одну статистику. По сгенерированным подвыборкам мы можем посчитать $B$ статистик и увидеть как наша неизвестная статистика распределена. Скорее всего, вам кажется это сложным. Давайте попробуем посмотреть на конкретный пример. Возможно, станет легче. Но это неточно.
Verizon — основная региональная телекоммуникационная компания (Incumbent Local Exchange Carrier, ILEC) в западной части США. В связи с этим данная компания обязана предоставлять сервис ремонта телекоммуникационного оборудования не только для своих клиентов, но и для клиентов других локальных телекоммуникационых компаний (Competing Local Exchange Carriers, CLEC). При этом в случаях, когда время ремонта оборудования для клиентов других компаний существенно выше, чем для собственных, Verizon может быть оштрафована.
Проверим правда ли время на ремонт своего оборудования существенно ниже чем время на ремонт других компаний.
library('dplyr')
library('rio') # тут import, который сам догадается до разделителей и даже excel прочтёт
df = import('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/verizon.csv')
# ILEC - внутренний клиент
# CLEC - внешние клиенты
head(df)
Time | Group |
---|---|
17.50 | ILEC |
2.40 | ILEC |
0.00 | ILEC |
0.65 | ILEC |
22.23 | ILEC |
1.20 | ILEC |
df %>% group_by(Group) %>% summarise(count = n(), mean = mean(Time), median = median(Time))
Group | count | mean | median |
---|---|---|---|
CLEC | 23 | 16.509130 | 14.33 |
ILEC | 1664 | 8.411611 | 3.59 |
qplot(df$Time)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Из-за этого средние будут не очень адекватной оценкой для типичного представителя выборки. Используя их, мы будем искуственно находиться на стороне компании. Поэтому гораздо уместнее сформулировать в терминах медиан, так как они нечувствительны к выбросам в отличие от средних.
$$ \begin{aligned} &H_0: \hspace{2mm} med_1 = med_2 \\ &H_1: \hspace{2mm} med_1 > med_2 \end{aligned} $$Беда в том, что для медиан как раз у нас нет хорошей статистики. Это весомый повод прибегнуть к бустрапу!
df %>% sample_n(size = 3, replace=TRUE)
Time | Group |
---|---|
9.48 | ILEC |
17.33 | ILEC |
0.25 | ILEC |
Сделали подвыборку размера $3$ с повторениями, потом по ней можно что-нибудь посчитать. Например, медианы. Закодим бустрап в виде двух функций. Первая будет делать B_sample
выборок.
get_bootstrap_samples <- function(data, B_sample){
samples = list() # лист, куда мы будем записывать выборки
for(i in 1:B_sample){
# генерируем выборки с повторениями соответствующую по размерам оригинальной
# и записываем её в лист
samples[[i]] = sample(data, size = length(data), replace=TRUE)
}
return(samples)
}
# неколько сбустрапированных выборк
# get_bootstrap_samples(df$Time, 3)
Вторая функция будет искать квантили сгерированого распределения.
stat_intervals <- function(stat, alpha){
return(quantile(stat, c(alpha/2,(1-alpha/2)), name = FALSE))
}
Отделим данные по внутренним и внешним клиентам друг от друга и построим доверительные интервалы для каждой из медиан. Будем бустрапировать $1000$ выборок.
ilec = df %>% filter(Group == 'ILEC')
clec = df %>% filter(Group == 'CLEC')
# Делаем 1000 выборок для первой таблички и считаем по каждой медиану
# функция lapply применяет функцию median к каждой выборке
# команда as.numeric сделает листы векторами
clec_median_scores = as.numeric(lapply(get_bootstrap_samples(clec$Time, 10000), median))
ilec_median_scores = as.numeric(lapply(get_bootstrap_samples(ilec$Time, 10000), median))
Смотрим на итоговые доверительный интервалы.
cat("95% confidence interval for the ILEC median repair time:", stat_intervals(ilec_median_scores, 0.05),'\n')
cat("95% confidence interval for the CLEC median repair time:", stat_intervals(clec_median_scores, 0.05))
95% confidence interval for the ILEC median repair time: 3.22 3.84 95% confidence interval for the CLEC median repair time: 5.45 20
По аналогии найдём доверительный интервал для разности медиан.
delta_median_scores = clec_median_scores - ilec_median_scores
cat("95% confidence interval for the diff median repair time:", stat_intervals(delta_median_scores, 0.05))
95% confidence interval for the diff median repair time: 2.02 16.66
Оценим вероятность того, что разность больше нуля.
sum(delta_median_scores > 0)/length(delta_median_scores)
Кажется, команию пришло время оштрафовать. Ещё один пример на бустрап мы посмотрим в тетрадке с поющими котами ближе к концу пары.
№
В 2012 году ряд авторов получил шнобелевскую премию по нейробиологии
Нужно было протестировать аппарат МРТ, обычно для этого берут шарик с маслом и сканируют его
Скучно, купили мёртвого лосося
Стали показывать ему фотки людей
Задача показать, что в голове лосося нет мозговой активности
Аппарат МРТ возвращает кучу данных. Один объект - воксель. Для исследования активности мозга в целом, надо провести множественное тестирование гипотез относительно каждого вокселя.
Оказалось, что у лосося есть рекация в мозгу на людей, написали статью
Уровень значимости: $P(\text{отвергнуть} H_0 \mid H_0 \text{верна}) = \alpha$
Выборка: $X_1, \ldots, X_n \sim iid \hspace{1mm} F_X(x, \theta)$
Нулевая гипотеза: $H_0 : \theta = \theta_0 $
Альтернатива: $H_1 : \theta \ne \theta_0 $
Статистика: $T(X_1, \ldots, X_n)$
Если $p$-значение $< \alpha$, гипотеза отвергается
Гипотеза сразу о двух параметрах
Всё усложняется: можно ошибиться в обеих гипотезах, ошибиться только в первой или только во второй
Накопление вероятности ошибки первого рода
Будем проверять каждую из двух гипотез на уровне значимости $\frac{\alpha}{2}$
Если гипотез $K$, то каждую проверяем на уровне $\frac{\alpha}{K}$
Такой подход называется поправкой Бонферрони
Если ввести поправку Бонферрони, мёртвый лосось остаётся мёртвым
Реальный вклад в науку, анализ работ по нйробиологии до 2010 года показал, что окло 40% не используют поправку при множественном тестировании гипотез
После публикации количество статей, где его не используют упало до 10%
Говорят, Джордж Р.Р. Мартин, автор цикла "Песнь Льда и Пламени", истребляет Старков: чаще убивает персонажей, относящихся к этому дому, чем персонажей других домов. В таблице ниже приведено количество персонажей, относящихся к тому или иному дому, упомянутых за первые $4$ книги, а так же количество погибших персонажей.
Дом | Упомянутые персонажи | Погибшие персонажи |
---|---|---|
House Stark | 72 | 18 |
House Lannister | 49 | 11 |
House Greyjoy | 41 | 12 |
Night's Watch | 105 | 41 |
Данные взяты из датасета по ссылке. Нужно проверить гипотезу об этом на $5\%$ уровне значимости.
Нам нужно проверить гипотезу о том, что доля смертей среди старков совпадает со всеми остальными долями против альтернативы, что хотя бы в одной из ситуаций смертность в доме Старков больше:
\begin{equation*} \begin{aligned} H_0: p_1 = p_2 \\ H_1: p_1 > p_2 \end{aligned} \end{equation*}Протестируем все гипотезы без коррекции на множественное тестирование, а потом введём её.
n_st = 72
m_st = 18
n_lan = 49
m_lan = 11
n_gr = 41
m_gr = 12
n_nw = 105
m_nw = 41
alpha = 0.05
Закодим тест для проверки гипотезы о равенстве долей в виде двух функций. Одна будет считать значение статистики, вторая находить p_value
.
proportions_diff_zstat = function(n1, m1, n2, m2){
p1 = m1/n1
p2 = m2/n2
P = (p1*n1 + p2*n2)/(n1 + n2)
z_stat = (p1 - p2)/sqrt(P * (1-P) * (1/n1 + 1/n2))
return(z_stat)
}
proportions_diff_ztest = function(z_stat, alternative = 'two-sided'){
if(alternative == 'two-sided'){
return(2 * (1 - pnorm(np.abs(z_stat))))
}
if(alternative == 'less'){
return(pnorm(z_stat))
}
if(alternative == 'greater'){
return(1 - pnorm(z_stat))
}
}
Проверяем три гипотезы!
# Для всех трёх случаев гипотеза о рпвенстве долей не отвергается.
proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_gr, m_gr),'greater')
proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_lan, m_lan),'greater')
proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_nw, m_nw),'greater')
pvalues = rep(0,3)
pvalues[1] = proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_gr, m_gr),'greater')
pvalues[2] = proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_lan, m_lan),'greater')
pvalues[3] = proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_nw, m_nw),'greater')
pvalues > 0.05/3 # не отвергается
Чтобы ввести коррекцию Бонферони, будем сравнивать с $\frac{\alpha}{3}$.
В чём минус такой поправки? Выше мы выписали неравенство
$$ P(\text{ошибочно отвергнуть } H_0 ) = P(\text{ошибочно отвергнуть хотя бы одну из гипотез}) \le \\ \le P(\text{ошибка в первой}) + P(\text{ошибка во второй}) = \alpha + \alpha = 2\alpha. $$В нём фигурировало две гипотезы. Если бы мы захотели бы по-честному выписать равенство, нужно было бы по формуле включений-исключений отнять вероятность пересечения событий. Если гипотез нужно проверить очень много, формула включе-ний исключений получается огромной. Мы в нашем неравенстве пользуемся только первым её членом. Получается, что реальный уровень значимости, на котором мы проверяем гипотезу оказывается намного меньше, чем заявленный нами $\alpha$. В идеале нам бы хотелось, чтобы вероятность совершить хотя бы одну ошибку первого рода была в точности равна $\alpha$.
Ещё раз, ещё раз: При использовании метода Бонферрони эта вероятность ограничивается гораздо более низкой величиной, чем$\alpha$. Это плохо, потому что перестраховываясь в отношении ошибки первого рода, мы неизбежно совершаем больше ошибок второго рода, то есть мощность такой статистической процедуры снижается.
Поправку Бонферрони очень просто применять, но у такой процедуры низкая мощность
Поправка Бонферрони берёт $\alpha_i$ одинаковыми
Можно улучшить процедуру, если брать $\alpha_i$ разными
Если $p_{(1)} \ge \alpha_1$, не отвергаем все нулевые гипотезы, иначе отвергаем $H_{(1)}$ и продолжаем
Если $p_{(2)} \ge \alpha_2$, не отвергаем все нулевые гипотезы, иначе отвергаем $H_{(2)}$ и продолжаем
Делаем так, пока не закончатся гипотезы
Нисходящая процедура, так как движемся по убыванию значимости
Если взять уровни значимости
$$\alpha_1 = \frac{\alpha}{m}, \alpha_2 = \frac{\alpha}{m-1}, \ldots, \alpha_i = \frac{\alpha}{m-i+1}, \ldots, \alpha_m = \alpha$$мы получим метод холма.
sort(pvalues)
seq(3,1,by=-1)
alpha/seq(3,1,by=-1)
sort(pvalues) > alpha/seq(3,1,by=-1) # гипотеза не отвергается
В поправках выше мы контролировали вероятность групповой ошибки, то есть вероятность совершить хотя бы одну ошибку первого рода
Когда проверяется слишком большое количество гипотез, можно допустить какое-то число ошибок первого рода ради повышения мощности процедуры, это позволит отвергнуть больше неверных гипотез
$\mbox{ }$ | число верных $H_i$ | число неверных $H_i$ |
---|---|---|
число не отвергнутых $H_i$ | U | T |
число отвернутых $H_i$ | V | S |
Если $p_{(m)} \le \alpha_m$, тогда все нулевые гипотезы отвергаются, иначе $H_m$ не отвергается и процедура продолжается
Если $p_{(m-1)} \le \alpha_{m-1}$, тогда все оставшиеся гипотезы отвергаются, иначе $H_{m-1}$ тоже не отвергаетмя и мы продолжаем процедуру
Делаем так, пока не закончатся гипотезы
Восходящая процедура, так как движемся по возростанию значимости, она помогает отвергать больше гипотез
Если для одних и тех же значений $\alpha_i$ построить восходящую и нисходящую процедуры, то восходящая процедура всегда будет отвергать не меньше гипотез, чем нисходящая.
Если взять уровни значимости
$$\alpha_1 = \frac{\alpha}{m}, \ldots, \alpha_i = \frac{i \cdot \alpha}{m}, \ldots, \alpha_m = \alpha$$мы получим метод Бенджамини-Хохберга. Тут в отличие от метода Холма уровни значимости возрастают по Гиперболе, а не линейно.
В R есть встроенная функция, позвоялющая сделать коррекцию на множественное тестирование гипотез
p.adjust(c(0.01,0.02,0.07), "bonferroni")
p.adjust(c(0.01,0.02,0.07), "holm")
p.adjust(c(0.01,0.02,0.03), "hochberg")