class: middle, left, inverse, title-slide .title[ # Тестирование статистических гипотез ] .subtitle[ ## Линейные модели… ] .author[ ### Марина Варфоломеева, Юта Тамберг, Вадим Хайтов ] .date[ ### Осень 2023 ] --- ## Тестирование гипотез - Выборочное распределение среднего значения - Как устроено тестирование гипотез - t-статистика - Применение t-теста - Статистическая значимость --- class: middle, center, inverse # Оценка параметров распределения по выборке --- ## Задание - Загрузите объект `population`. Это вектор, характеризующий размеры особей, в некоторой гипотетической популяции. ```r load("data/population.RData") ``` - Сделайте из этой совокупности случайную выборку в 5 измерений. ```r my_sample <- sample(x = population, size = 5) my_sample ``` ``` [1] 175.4 167.9 173.7 179.4 172.8 ``` - Определите среднее значение и среднеквадратичное отклонение. ```r mean(my_sample) ``` ``` [1] 173.8 ``` ```r sd(my_sample) ``` ``` [1] 4.175 ``` - Проведите такую выборку несколько раз. --- ## Выборочные оценки - Изучая интересующую нас систему, мы не знаем значений параметров распределения, которое описывает поведение вероятности в **генральной совокупности** (*population*). - Мы проводим **оценку** (*estimate*) этих параметров, используя **выборку** (*sample*). - Выборка имеет ограниченный **объем**(*N*, или *n*). - Выборка попадает в анализ в виде вектора значений `\(\{x_1, x_2, x_3, \dots , x_n\}\)` --- ## Выборочные оценки Выборка может быть описана с помощью "5 number statistic": -- - Min - Q1 - Q2 (Median) - Q3 - Max --- ## Выборочные оценки Другой способ описания выборки -- - Среднее значение: `\(\bar{x}=\frac{\Sigma{x_i} }{n}\)` -- - Стандартное отклонение: `\(sd = \sqrt{\frac{\Sigma{(x_i - \bar{x})^2}}{n-1} }\)` -- Если `\(x_i \in N(\mu, \sigma)\)`, то -- $$ \bar{x} \xrightarrow[n\to\infty]{} \mu $$ -- $$ sd \xrightarrow[n\to\infty]{} \sigma $$ --- class: middle, center, inverse # Выборочная оценка среднего ??? Что мы можем сказать о среднем значении в генеральной совокупности, если у нас есть всего одна выборка? Центральная предельная теорема говорит, что если мы возьмем достаточно большую выборку из генеральной совокупности, то среднее значение будет нормально распределено. Особенно важно, что это правда даже если признак в совокупности имеет другое распределение. Повторные выборки. Зависимость точности оценки от объема выборки. Ошибка среднего SE (это иллюстрация теоремы центрального предела). --- ## Как можно судить о свойствах генеральной совокупности по выборке? __Центральная предельная теорема__ (ЦПТ) говорит, что если мы возьмем достаточно большую выборку из генеральной совокупности, то среднее значение будет нормально распределено с параметрами `\(\mu_{\bar x}\)` и `\(\sigma _{\bar{x}}\)`: `$$\bar X \sim N (\mu_{\bar x}, \sigma_{\bar x})$$` При чем `\(\sigma_{\bar x} = \frac{\sigma}{\sqrt{n}}\)`. <br/> __Важно__: это будет так при *больших объемах выборки* ( `\(n > 30\)`, или даже `\(n > 100\)`), даже если `\(x\)` в генеральной совокупности не подчиняется нормальному распределению. <br/> Давайте проверим на опыте, так ли это. --- ## Цена алмазов Представим, что данные об алмазах из датасета `diamonds` (пакет `gplot2`) — это генеральная совокупность. Перед вами распределение цены алмазов. Давайте будем брать из этого распределения выборки и оценивать по ним среднее значение. .pull-left[ <!-- --> ] .pull-right[  .tiny[Diamond Age by Steve Jurvetson on Flickr <!-- https://flic.kr/p/eRNcR -->] ] --- .pull-left[ <!-- --> ] .pull-right[ ## Средние в выборках Средние в выборках отличаются от среднего в генеральной совокупности. Если взять много выборок определенного размера, можно построить распределение выборочных средних. <br/> Как изменится форма распределения выборочных средних при изменении объема выборки? ] --- .pull-left[ <!-- --> ] .pull-right[ ### Распределение выборочных средних при разных объемах выборки `$$\bar X \sim N (\mu_{\bar x}, \sigma_{\bar x})$$` `\(\mu_{\bar x} = \mu\)` — среднее значение выборочных средних стремится к среднему в генеральной совокупности. `\(\sigma_{\bar x} = \sigma / \sqrt{n}\)` — стандартное отклонение в `\(\sqrt{n}\)` раз меньше стандартного отклонения в генеральной совокупности. `\(\sigma_{\bar x}\)` называют стандартной ошибкой среднего и обозначают `\(SE _{\bar{x}}\)`. ] --- ## Центральная предельная теорема очень важна в статистике .pull-left[ <!-- --> ] .pull-right[ `$$\bar X \sim N (\mu, \sigma / \sqrt{n})$$` Пользуясь ее выводами, мы сможем: - строить доверительные интервалы - тестировать гипотезы ] --- class: middle, center, inverse # Доверительный интервал --- ## Если выполняется центральная предельная теорема... <br/><br/> .pull-left-60[ <!-- --> ] .pull-right-40[ Было `$$\bar X \sim N(\mu, \sigma/ \sqrt{n})$$` После стандартизации: `$$\frac{\bar X - \mu}{\sigma / \sqrt{n}} \sim N(0, 1)$$` Стандартизованное распределение выборочных средних — это стандартное нормальное распределение. ] --- ## Доверительный интервал <br/>из нормального распределения .pull-left-60[ <!-- --> ] .pull-right-40[ --- это интервал, в который попадает заданный процент выборочных средних. В 95% доверительный интервал попадает выборочное среднее в 95% *повторных* выборок. <br/> Как найти этот интервал? ] --- ## Доверительный интервал <br/>из нормального распределения .pull-left-60[ <!-- --> ] .pull-right-40[ `$$\bar {x} \pm z_{\alpha} \cdot \sigma / \sqrt{n}$$` Чтобы найти границы 95% доверительного интервала, нужно найти квантили стандартного нормального распределения, которые соответствуют вероятностям 0.025 и 0.975 ```r qnorm(p = c(0.025, 0.975)) ``` ``` [1] -1.96 1.96 ``` `\(z_{0.05} = 1.96\)` 95% выборочных средних в повторных выборках будут лежать в пределах `\(\pm 1.96\)` стандартных ошибок вокруг среднего значения. ] --- ## Условия применимости нормального распределения для доверительного интервала 1.Должна быть известна `\(\sigma\)` в генеральной совокупности. 2.Должны выполняться условия, при которых справедлива ЦПТ: - Наблюдения в выборке должны быть независимы друг от друга. - Большой объем выборки **или** нормальное распределение `\(x\)` --- ## Если `\(\pmb \sigma\)` не известна Если `\(\sigma\)` в генеральной совокупности не известна, ее можно оценить по выборочному стандартному отклонению `\(sd\)`. `$$\sigma / \sqrt{n} \approx sd/\sqrt{n}$$` После стандартизации: `$$\frac{\bar X - \mu}{SE_{\bar x}} = \frac{\bar X - \mu}{sd / \sqrt{n}} \sim t_{df = n - 1}$$` стандартизованное распределение выборочных средних подчиняется `\(t\)`-распределению с числом степеней свободы `\(df = n - 1\)` --- ## _t_-распределение, или распределение Стьюдента .pull-left-55[ <!-- --> - Симметричное колоколообразное распределение с толстыми хвостами. - Единственный параметр — число степеней свободы (для доверительного интервала `\(df = n - 1\)`). - При увеличении объема выборки `\(t\)`-распределение приближается к нормальному. ] .pull-right-45[  William Sealy Gosset ] --- ## Доверительный интервал из _t_-распределения .pull-left-55[ <!-- --> ] .pull-right-45[ Обязательно используется, если: - Объем выборки мал. - `\(\sigma\)` не известна. `$$\bar {x} \pm t_{\alpha, df} \cdot s / \sqrt{n}$$` `\(df = n - 1\)` <br/><br/><br/><br/><br/><br/> ] Условия применимости Выполняются условия, при которых справедлива ЦПТ: - Наблюдения в выборке независимы друг от друга. - Большой объем выборки и нет "выбросов" **или** нормальное распределение `\(x\)` --- ## Смысл 95% доверительного интервала .pull-left[ <!-- --> ] .pull-right[ - Для примера возьмем данные из объекта `population`. Здесь известно, что `\(\mu = 170\)`. - Симулируем взятие многих выборок из этой генеральной совокупности. - Среднее в генеральной совокупности — это фиксированная величина (она либо попала в интервал, либо нет. - Доверительный интервал — случайная величина. - В *повторных выборках* одинакового объема `\(\approx 95\%\)` всех доверительных интервалов "накроют" истинное среднее значение. ] --- ## Расчет и изображение доверительного интервала в R ```r ## library(ggplot2) data("diamonds") # цена бриллиантов хорошего качества огранки good <- diamonds$price[diamonds$cut == "Good"] .mean <- mean(good) # выборочное среднее .n <- length(good) # объем выборки SE <- sd(good)/ sqrt(.n) # стандартная ошибка t_crit <- qt(p = 0.975, df = .n - 1) # критич. зн. t для данного n и p = 0.95 err <- t_crit * SE # предел погрешности err ``` ``` [1] 103 ``` ```r # Границы доверительного интервала .mean - err ``` ``` [1] 3826 ``` ```r .mean + err ``` ``` [1] 4032 ``` Можем записать среднюю цену бриллиантов хорошей огранки и ее доверительный интервал: `\(3928.9 \pm 103 (\pm SE)\)` --- ## Строим доверительные интервалы в ggplot ```r ## theme_set(theme_bw()) ggplot(data = diamonds, aes(x = cut, y = price)) + stat_summary(geom = 'pointrange', fun.data = mean_cl_normal) ``` <!-- --> --- ## Задание 1 Посчитайте среднюю цену и доверительный интервал для бриллиантов с такими свойствами (*одновременно*): - качество огранки (`cut`) идеальное (`Ideal`) - прозрачность (`clarity`) наивысшая (`IF`) Постройте график средней цены с доверительными интервалами для бриллиантов разного качества огранки и прозрачности. --- ## Решение Вычислим границы доверительного интервала ```r ideal <- diamonds$price[diamonds$cut == 'Ideal' & diamonds$clarity == 'IF'] .mean <- mean(ideal) # выборочное среднее .n <- length(ideal) # объем выборки, если нет NA SE <- sd(ideal)/ sqrt(.n) # стандартная ошибка t_crit <- qt(p = 0.975, df = .n - 1) # критич. зн. t для данного n и p = 0.95 err <- t_crit * SE # предел погрешности err ``` ``` [1] 187.1 ``` ```r # Границы доверительного интервала .mean - err ``` ``` [1] 2086 ``` ```r .mean + err ``` ``` [1] 2460 ``` --- ## Решение: cпособ 1 (некрасивый, требует доработки) Теперь можно построить график средней цены с доверительными интервалами для бриллиантов разного качества огранки и прозрачности. Для этого нужно к предыдущему графику добавить информацию о прозрачности. ```r ggplot(data = diamonds, aes(x = cut, y = price, colour = clarity)) + stat_summary(geom = 'pointrange', fun.data = mean_cl_normal) ``` <!-- --> --- ## Решение: cпособ 1 (доработка) Теперь можно построить график средней цены с доверительными интервалами для бриллиантов разного качества огранки и прозрачности. Для этого нужно к предыдущему графику добавить информацию о прозрачности. ```r ggplot(data = diamonds, aes(x = cut, y = price, colour = clarity)) + stat_summary(geom = 'pointrange', fun.data = mean_cl_normal, position = position_dodge(width = 0.5)) ``` <!-- --> --- ## Решение: cпособ 2 (удовлетворительный) ```r ggplot(data = diamonds, aes(x = cut, y = price, colour = clarity)) + stat_summary(geom = 'pointrange', fun.data = mean_cl_normal) + facet_wrap(~ clarity) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) ``` <!-- --> Что странного в этой картинке? -- Осторожно! Группы могут быть неоднородны: в каждой группе варьирует размер алмазов. --- ## Статистика по группам при помощи пакета `dplyr` На примере доверительных интервалов: ```r library(dplyr) smr_diamonds <- diamonds %>% # данные group_by(cut, clarity) %>% # группируем summarize( .mean = mean(price), # выборочное среднее .n = n(), # объем выборки SE = sd(price) / sqrt(.n), # стандартная ошибка t_crit = qt(p = 0.975, df = .n - 1), # критич. зн. t err = t_crit * SE, # предел погрешности lower = .mean - err, # нижняя граница доверительного интервала upper = .mean + err # верхняя граница доверительного интервала ) head(smr_diamonds) ``` ``` # A tibble: 6 × 9 # Groups: cut [1] cut clarity .mean .n SE t_crit err lower upper <ord> <ord> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> 1 Fair I1 3704. 210 214. 1.97 422. 3282. 4125. 2 Fair SI2 5174. 466 182. 1.97 358. 4816. 5531. 3 Fair SI1 4208. 408 162. 1.97 319. 3889. 4527. 4 Fair VS2 4175. 261 219. 1.97 431. 3744. 4605. 5 Fair VS1 4165. 170 276. 1.97 546. 3619. 4711. 6 Fair VVS2 3350. 69 364. 2.00 727. 2623. 4077. ``` --- class: middle, center, inverse # Как устроено тестирование гипотез --- ## Сравнение выборок Различия между выборками не всегда видны невооружённым глазом.   .tiny[tres caracoles by Alberto Villen on Freeimages.com] --- ## Нулевая и альтернативная гипотезы Это первый шаг   .tiny[tres caracoles by Alberto Villen on Freeimages.com] -- - Нулевая гипотеза `\(H_0\)` чаще всего формулируется как **отсутствие различий** между сравниваемыми объектами. Например: Улитки из обеих популяций одинакового размера - Альтернативная гипотеза `\(H_A\)` формулируется как **присутствие различий**, она обратна нулевой гипотезе, т.е. включает все остальные случаи. Например: Улитки из обеих популяций разного размера. --- ## Нулевая и альтернативная гипотезы — это "два мира" Вне зависимости от нас, реальность может находиться в одном из двух состояний: .pull-left[ - `\(H_0\)` верна, улитки одинаковы   ] .pull-right[ - `\(H_0\)` неверна, улитки различаются   ] <br /> -- После статистического теста мы принимаем решение о том, принять или отвергнуть `\(H_0\)`. Но это решение не обязательно окажется верным. Возможно четыре исхода: .pull-left[ В мире где улитки одинаковы <br/>( `\(H_0\)` верна) мы можем: - принять `\(H_0\)` (верное решение), - отвергнуть `\(H_0\)` (ошибка). ] .pull-right[ В мире где улитки различаются <br/>( `\(H_A\)` верна), мы можем: - принять `\(H_0\)` (ошибка), - либо отвергнуть `\(H_0\)` (верное решение). ] --- ## Верные и неверные решения .pull-left[ **Ошибка I рода: нашли то, чего нет** ] .pull-right[ **Ошибка II рода: не нашли то, что было** ] | | `\(H_0\)` верна | `\(H_0\)` неверна | |:-----:|:-----:|:-----:| | Отклонить H0 | Ошибка I рода с вероятностью <span class="orange">&alpha;</span></br>Ложно-положительный результат | Верно | | Сохранить H0 | Верно | Ошибка II рода с вероятностью <span class= "blue">&beta;</span> </br> Ложно-отрицательный результат | -- Пока мы будем говорить о том, как контролируют уровень **Ошибок I рода**. Про **Ошибки II рода** и **Power analysis** читайте дополнительную литературу. --- ## Тестирование гипотез: Тестовые статистики -- #### 1. Формулируем нулевую и альтернативную гипотезы. Гипотезы выражаются математически в виде тестовых статистик. На этом этапе мы делаем определенные допущения. -- #### 2. Проверяем __условия применимости__ тестовой статистики. -- #### 3. По реальным данным вычисляем __эмпирическое значение тестовой статистики__. -- Дальше мы должны ответить на вопрос: **Насколько вероятно получить _такое или более экстремальное_ эмпирическое значение, если верна нулевая гипотеза `\(H_0\)`?** #### 4. Строим теоретическое распределение тестовой статистики для случая, когда верна `\(H_0\)`, и оцениваем по нему уровень значимости. -- #### 5. Решаем сохранить или отвергнуть `\(H_0\)`. Увы, мы не сможем узнать, какая гипотеза верна, но поймем, насколько с ней согласуются исходные данные. --- class: middle, center, inverse # Одновыборочный *t*-тест --- ## Размер кладки черепах Разберемся с одновыборочным `\(t\)`-тестом на вымышленном примере. Представьте, что в [одной статье](https://edis.ifas.ufl.edu/publication/UW441#:~:text=Reproductive%20rate%3A%20Clutch%20sizes%20range,to%20100%20years%20in%20captivity.) сказано, что средняя плодовитость черепах определенного вида — 8 яиц в кладке. В вашей выборке из `\(35\)` черепах — `\(\bar x = 9\)`, `\(sd = 2\)`. Отличается ли приведенное в статье значение от того, что наблюдается в изученной популяции. Можно ли считать, что авторы статьи изучали что-то иное (другой вид, или популяцию из иных условий). <!-- Отличается ли реальная плодовитость в обследованной вами популяции черепах от того, что указано в статье? -->  <small>Gopher Tortoise by Judy Gallagher on Flickr</small> <!-- https://flic.kr/p/Q2ZozS --> --- ## Одновыборочный t-тест - `\(H_0: \mu = \mu_0\)` — Реальная средняя плодовитость черепах такая, как в статье. - `\(H_A: \mu \ne \mu_0\)` — Средняя плодовитость отличается от того, что написано в статье. `\(\mu_0\)` — это какое-то конкретное значение. В нашей задаче это — 8 яиц в кладке. <br/> `$$t = \cfrac{\bar x - \mu}{ sd / \sqrt{n} }$$` Если выполняется ЦПТ, то одновыборочная `\(t\)`-статистика подчиняется `\(t\)`-распределению с числом степеней свободы `\(df = n - 1\)`. Условия применимости: - Наблюдения в выборке должны быть независимы друг от друга. - Объем выборки достаточно велик **или** величины нормально распределены. ### Задание 2 Проверьте условия применимости t-теста. Вычислите t и p. --- ## Проверяем, нормально ли распределение <!-- --> -- Нет оснований для сомнения, что наша выборка взята из генеральной совокупности с распределением, которое отличается от нормального. --- ## Вычислим наблюдаемое значение *t*-статистики `$$t = \cfrac{\bar x - \mu}{ sd / \sqrt{n} }$$` Средняя плодовитость в выборке из 35 черепах `\(\bar x = 9\)`, стандартное отклонение `\(s = 2\)`. В статье указана плодовитость `\(8\)`. `$$t = \cfrac{9 - 8}{ 2 / \sqrt{35} } = 2.96$$` --- ## Насколько это значение согласуется с `\(H_0\)`? .pull-left[ <!-- --> ] .pull-right[ При `\(H_0\)` значение `\(t\)` будет близко к нулю. Насколько необычны значения t меньше или больше `\(\pm 2.96\)`? ] --- ## Уровень значимости (_p_-value) .pull-left[ <!-- --> ] .pull-right[ __Уровень значимости__ — это вероятность получить значение `\(t\)` меньше или больше данного, если бы `\(H_0\)` была справедлива. ] --- ## Уровень значимости (_p_-value) .pull-left[ <!-- --> ] .pull-right[ Можно вычислить значение `\(p\)` ```r 2 * ( 1 - pt(2.96, df = 35 - 1) ) ``` ``` [1] 0.005571 ``` Если бы плодовитость не отличалось от указанной в статье, получить `\(t\)` меньше или больше `\(2.96\)` можно было бы с вероятностью `\(p = 0.006\)`. ] --- ## Критический уровень значимости .pull-left[ <!-- --> ] .pull-right[ Критический уровень значимости `\(\alpha\)` — это порог для принятия решений. Обычно используют `\(\alpha = 0.05\)`. `\(\alpha\)` - это допустимая вероятность **Ошибки I рода** (найти отличия там, где их нет). Если `\(p \le \alpha\)` — отвергаем `\(H_0\)` и принимаем `\(H_A\)`. Если `\(p > \alpha\)` — сохраняем `\(H_0\)`, не можем ее отвергнуть. ] --- ## Принимаем решение .pull-left[ <!-- --> ] .pull-right[ Мы получили `\(p < \alpha\)`, поэтому отвергаем `\(H_0\)` и принимаем `\(H_A\)`. Указанная в статье плодовитость **статистически значимо** отличается от плодовитости в изученной нами популяции. Мы изучали что-то иное (другую генеральную совокупность)! -- *At!* Лучше не использовать популярное словосочетание "*достоверно отличается*". Достоверные события - это события, вероятность которых равна 1. ] --- ## Заблуждения о _p_-values -- Правда ли, что `\(p\)` — вероятность того, что верна сама `\(H_0\)`? -- - Нет! Значение `\(p\)` всегда считается __при условии, что `\(H_0\)` верна__. Но `\(H_0\)` может быть неверна. <br/> -- Правда ли, что `\(p\)` — это вероятность получить такое значение статистики при справедливой `\(H_0\)`? -- - Нет! Вероятность вычисляется как площадь под участком кривой. Конкретное значение статистики — это точка и под ней нет площади. <br/> -- Правда ли, что если `\(p > 0.05\)`, то различий между группами на самом деле нет? -- - Нет! Это значит, что имеющиеся данные не позволяют отвергнуть нулевую гипотезу. Различия могут быть. -- *p*-value - это всего лишь вероятность того, что *в мире, где справедлива `\(H_0\)`,* при повторных выборках из той же генеральной совокупности значение выбранной статистики (например *t*) будет больше или равно наблюдаемого в нашей выборке. --- class: middle, center, inverse # Двухвыборочный t-тест --- ## Гипотезы в двухвыборочном `\(t\)`-тесте и тестовая статистика `\(H_0: \mu_1 - \mu_2 = 0\)` — средние значения не различаются в двух группах `\(H_A: \mu_1 - \mu_2 \ne 0\)` — средние значения различаются <br/> Т.е. нас интересует __разность выборочных средних__. Ее ожидаемое значение при `\(H_0\)` будет 0. t-тест в общем виде выглядит так `$$t=\frac{\text{Наблюдаемая величина - Ожидаемое значение}}{\text{Стандартная ошибка}}$$` --- ## t-тест и его разновидности Двухвыборочный t-тест используется для проверки значимости различий между средними `$$t\;= \;\frac{\bar{x}_1 - \bar{x}_2}{SE_{\bar{x}_1 - \bar{x}_2}}$$` -- `\(SE_{\bar{x}_1 - \bar{x}_2}\)` — стандартная ошибка разности двух средних, может рассчитываться по-разному - t-тест Стьюдента — если считать, что дисперсии в группах равны - t-тест Уэлча — если считать, что дисперсии могут быть разными --- ## Стандартная ошибка разности средних в t-тесте Стьюдента Student 1908 Если группы независимы и дисперсии в них равны, то по центральной предельной теореме `$$SE_{\bar{x}_1 - \bar{x}_2} = \sqrt{ \frac{\sigma^2}{n_{1}} + \frac{\sigma^2}{n_{2}}} \approx \sqrt{ \frac{sd^2}{n_{1}} + \frac{sd^2}{n_{2}}}$$` `\(s^2 = \cfrac{(n_1 - 1)sd^2_1 + (n_2 - 1)sd^2_2}{n_1 + n_2 - 2}\)` — это __обобщенная дисперсия__ по двум выборкам. Результирующая `\(t\)`-статистика `$$t\;= \;\frac{\bar{x}_1 - \bar{x}_2}{\sqrt{SE_1^2 + SE_2^2}}$$` подчиняется `\(t\)`-распределению с `\(df = n_1 + n_2 - 2\)`. -- Осторожно, равенство дисперсий в группах — это часто нереалистичное предположение! --- ## Cтандартная ошибка разности средних в t-тесте Уэлча Если группы независимы и дисперсии в них неизвестны, то получается `$$SE_{\bar{x}_1 - \bar{x}_2} = \sqrt{ \frac{\sigma^2_{1}}{n_{1}} + \frac{\sigma^2_{2}}{n_{2}}} \approx \sqrt{ \frac{sd^2_{1}}{n_{1}} + \frac{sd^2_{2}}{n_{2}}}$$` -- Проблема в том, что эта величина __лишь приблизительно следует t-распределению__, если считать его степени свободы как обычно для двух групп `\(df = n_1 + n_2 - 2\)`. --- ## Приблизительное число степеней свободы Можно рассчитать по уравнению Уэлча-Саттеруэйта. Это решит проблему. `$$df_{ Welch-Satterthwaite} \approx \cfrac {\bigg(\cfrac{sd^2_{1}}{n_{1}} + \cfrac{sd^2_{2}}{n_{2}}\bigg)^2} {\cfrac{1}{n_{1} - 1}\bigg(\cfrac {sd_{1}^2} {n_{1}}\bigg)^2 + \cfrac{1}{n_{2} - 1}\bigg(\cfrac {sd_{2}^2} {n_{2}}\bigg)^2}$$` t-тестом Уэлча можно пользоваться, даже если дисперсии равны (это более общее решение). Он немного консервативнее, чем тест Стьюдента. --- ## Условия применимости двухвыборочного t-теста Почти такие же, как условия справедливости ЦПТ - Наблюдения независимы друг от друга. - Выборки независимы друг от друга (новое условие). - Объем выборки достаточно велик или величины нормально распределены. --- class: middle, center, inverse # Двухвыборочный t-тест в R --- ## Пример: Гормоны и артериальная гипертензия Синдром Кушинга — это нарушения уровня артериального давления, вызванные гиперсекрецией кортизола надпочечниками. В датасете `Cushings` (пакет `MASS`) записаны данные о секреции двух метаболитов при разных типах синдрома (данные из кн. Aitchison, Dunsmore, 1975). - `Tetrahydrocortisone` — секреция тетрагидрокортизона с мочой (мг/сут.) - `Pregnanetriol` — секреция прегнантриола с мочой (мг/сут.) - `Type` — тип синдрома: - `a` — аденома - `b` — двусторонняя гиперплазия - `c` — карцинома - `u` — не известно Давайте сравним секрецию тетрагидрокортизона при аденоме и двусторонней гиперплазии надпочечников. Различается ли она? --- ## Открываем данные ```r library(MASS) data("Cushings") ``` Все ли правильно открылось? ```r head(Cushings) ``` ``` Tetrahydrocortisone Pregnanetriol Type a1 3.1 11.70 a a2 3.0 1.30 a a3 1.9 0.10 a a4 3.8 0.04 a a5 4.1 1.10 a a6 1.9 0.40 a ``` ```r str(Cushings) ``` ``` 'data.frame': 27 obs. of 3 variables: $ Tetrahydrocortisone: num 3.1 3 1.9 3.8 4.1 1.9 8.3 3.8 3.9 7.8 ... $ Pregnanetriol : num 11.7 1.3 0.1 0.04 1.1 0.4 1 0.2 0.6 1.2 ... $ Type : Factor w/ 4 levels "a","b","c","u": 1 1 1 1 1 1 2 2 2 2 ... ``` --- ## Знакомимся с данными Есть ли пропущенные значения? ```r colSums(is.na(Cushings)) ``` ``` Tetrahydrocortisone Pregnanetriol Type 0 0 0 ``` <br/> Каковы объемы выборки в каждой группе? ```r table(Cushings$Type) ``` ``` a b c u 6 10 5 6 ``` Обратите внимание, объемы выборок **маленькие**. --- ## Проверяем условия применимости... 1.Наблюдения независимы друг от друга? - Да, независимы. Это случайная выборка. <br/> 2.Выборки независимы друг от друга? - Да, независимы. В группах разные люди (естественно, т.к. тип синдрома у человека может быть только какой-то один). <br/> 3.Объем выборки достаточно велик или величины нормально распределены? - Объем выборки мал ```r table(Cushings$Type) ``` ``` a b c u 6 10 5 6 ``` Нужно проверить форму распределения в обеих группах. --- ## Нормально ли распределены концентрации тетрагидрокортизона в группах? ```r library(car) qqPlot(Cushings$Tetrahydrocortisone[Cushings$Type == 'a']) qqPlot(Cushings$Tetrahydrocortisone[Cushings$Type == 'b']) ``` <!-- --> Удовлетворительно. При таких малых объемах выборки сложно ожидать лучшего. <!-- Будем считать, что можно аппроксимировать концентрацию тетрагидрокортизона нормальным распределением. --> --- ## Двухвыборочный t-тест в R: способ 1. Сравним секрецию тетрагидрокортизона при помощи **двухвыборочного** t-теста. ``` t.test(x = значения_в_гр.1, y = значения_в_гр.2) ``` ```r tt <- t.test(x = Cushings$Tetrahydrocortisone[Cushings$Type == 'a'], y = Cushings$Tetrahydrocortisone[Cushings$Type == 'b']) tt ``` ``` Welch Two Sample t-test data: Cushings$Tetrahydrocortisone[Cushings$Type == "a"] and Cushings$Tetrahydrocortisone[Cushings$Type == "b"] t = -4.1, df = 11, p-value = 0.002 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -7.988 -2.438 sample estimates: mean of x mean of y 2.967 8.180 ``` --- ## Опишем результаты ``` Welch Two Sample t-test data: Cushings$Tetrahydrocortisone[Cushings$Type == "a"] and Cushings$Tetrahydrocortisone[Cushings$Type == "b"] t = -4.1, df = 11, p-value = 0.002 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -7.988 -2.438 sample estimates: mean of x mean of y 2.967 8.180 ``` - Секреция тетрагидрокортизона значимо различается у пациентов с аденомой и двусторонней гиперплазией надпочечников ( `\(t_{10.69} = -4.15\)`, `\(p = <0.05\)`) <br/> Можно указать в скобках не сравнение с `\(\alpha\)`, а само значение `\(p\)`: ( `\(t_{10.69} = -4.15\)`, `\(p = 0.0017\)`). Только не надо безумствовать и указывать слишком много знаков... --- ## Задание 3: Двухвыборочный t-тест в R, способ 2. Перепишите вызов функции t.test с использованием другого шаблона вызова (с использованием формулы). `t.test(formula = что_зависит ~ группы, data = данные)` Второй вариант синтаксиса можно использовать только, если у вас есть только две группы. Если групп больше, то можно отфильтровать нужные с помощью логического вектора. ``` t.test(formula = что_зависит ~ группы, data = данные, subset = логический_вектор) ``` --- ## Решение: Двухвыборочный t-тест в R, способ 2. `t.test(formula = что_зависит ~ группы, data = данные)` Второй вариант синтаксиса можно использовать только, если у вас есть только две группы. Если групп больше, то можно отфильтровать нужные с помощью логического вектора. ``` t.test(formula = что_зависит ~ группы, data = данные, subset = логический_вектор) ``` ```r t.test(formula = Tetrahydrocortisone ~ Type, data = Cushings, subset = Cushings$Type %in% c('a', 'b')) ``` Результаты точно такие же, естественно. --- ## Задание 4 Посмотрите структуру результатов (`tt`) при помощи функции `str()` и извлеките из них: - степени свободы - уровень значимости - значение t-критерия --- ## Решение: Что спрятано в результатах? Как называются отдельные элементы результатов можно узнать посмотрев их структуру при помощи функции `str()` ```r str(tt) ``` ``` List of 10 $ statistic : Named num -4.15 ..- attr(*, "names")= chr "t" $ parameter : Named num 10.7 ..- attr(*, "names")= chr "df" $ p.value : num 0.00172 $ conf.int : num [1:2] -7.99 -2.44 ..- attr(*, "conf.level")= num 0.95 $ estimate : Named num [1:2] 2.97 8.18 ..- attr(*, "names")= chr [1:2] "mean of x" "mean of y" $ null.value : Named num 0 ..- attr(*, "names")= chr "difference in means" $ stderr : num 1.26 $ alternative: chr "two.sided" $ method : chr "Welch Two Sample t-test" $ data.name : chr "Cushings$Tetrahydrocortisone[Cushings$Type == \"a\"] and Cushings$Tetrahydrocortisone[Cushings$Type == \"b\"]" - attr(*, "class")= chr "htest" ``` --- ## Решение: Можно получить элементы результатов в виде отдельных чисел ```r tt$parameter # степени свободы ``` ``` df 10.69 ``` ```r tt$p.value # уровень значимости ``` ``` [1] 0.001719 ``` ```r tt$statistic # значение t-критерия ``` ``` t -4.15 ``` --- ## График со средними и доверительными интервалами ```r ggplot(data = Cushings[Cushings$Type %in% c('a', 'b'), ], aes(x = Type, y = Tetrahydrocortisone)) + stat_summary(fun.data = mean_cl_normal) ``` <!-- --> --- class: middle, center, inverse # Двухвыборочный t-тест (самостоятельная работа) --- ## Задание 5 Файл `aml.csv` содержит данные о влиянии регулярной химиотерапии на продолжительность ремиссии. Прочитаем эти данные ```r rem <- read.csv(file = "data/aml.csv", header = TRUE) head(rem) ``` ``` time group 1 9 1 2 13 1 3 13 1 4 18 1 5 23 1 6 28 1 ``` - В переменной `time` представлена продолжительность ремиссии в днях. - `group` указывает, к какой экспериментальной группе принадлежал пациент. В группе 1 проводилась регулярная химиотерапия, в группе 2 - нет. Ваша задача сравнить эти группы с помощью t-теста. --- ## Решение: Разведочный анализ ```r str(rem) ``` ``` 'data.frame': 23 obs. of 2 variables: $ time : int 9 13 13 18 23 28 31 34 45 48 ... $ group: int 1 1 1 1 1 1 1 1 1 1 ... ``` ```r # Делаем group фактором rem$group <- factor(rem$group, labels = c("therapy", "control")) # Пропущенные значения? colSums(is.na(rem)) ``` ``` time group 0 0 ``` ```r # Объемы выборок? table(rem$group) ``` ``` therapy control 11 12 ``` --- ## Решение: Проверка условий применимости ```r # Выборки малы, поэтому проверяем на нормальность и отсутствие выбросов op <- par(mfrow = c(1, 2)) qqPlot(rem$time[rem$group == "therapy"]) ``` ``` [1] 11 1 ``` ```r qqPlot(rem$time[rem$group == "control"]) ``` <!-- --> ``` [1] 12 11 ``` ```r par(op) ``` --- ## Решение: Проверка условий применимости Есть отскакивающее значение в группе `therapy`, в этой группе оно 11 по порядку (видно на `qqPlot()`). ```r rem$time[rem$group == "therapy"] ``` ``` [1] 9 13 13 18 23 28 31 34 45 48 161 ``` ```r rem$time[rem$group == "therapy"][11] ``` ``` [1] 161 ``` Пусть в этом случае это вполне реальное наблюдение, но это значение придется исключить, т.к. t-тест не устойчив к выбросам. --- ## Решение: Двухвыборочный t-тест ```r # удалено наблюдение со значением 161 tt <- t.test(formula = time ~ group, data = rem, subset = rem$time != 161) ``` Или то же самое: ```r # удалено 11 наблюдение в группе `therapy` tt <- t.test(x = rem$time[rem$group == 'therapy'][-11], y = rem$time[rem$group == 'control']) ``` --- ## Решение: Результаты двухвыборочного t-теста ```r tt ``` ``` Welch Two Sample t-test data: rem$time[rem$group == "therapy"][-11] and rem$time[rem$group == "control"] t = 0.83, df = 20, p-value = 0.4 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -7.486 17.386 sample estimates: mean of x mean of y 26.20 21.25 ``` - Не найдено значимых различий продолжительности ремиссии у пациентов после регулярной лучевой терапии и обычного лечения ( `\(t_{19.68} = 0.83\)`, `\(p = 0.416\)`). --- ## Решение: График со средними и доверительными интервалами ```r ggplot(data = rem[rem$time != 161, ], aes(x = factor(group), y = time)) + stat_summary(fun.data = mean_cl_normal) ``` <!-- --> --- class: middle, center, inverse # Парный t-тест (факультативно) --- ## Пример: Снотворное В датасете `sleep` содержатся данные об увеличении продолжительности сна по сравнению с контролем после применения двух снотворных препаратов (Cushny, Peebles, 1905, Student, 1908) Одинаково ли два снотворных влияют на увеличение продолжительности сна? ```r data(sleep) head(sleep) ``` ``` extra group ID 1 0.7 1 1 2 -1.6 1 2 3 -0.2 1 3 4 -1.2 1 4 5 -0.1 1 5 6 3.4 1 6 ``` ```r # str(sleep) ``` --- ## Проверяем условия применимости... 1.Наблюдения независимы друг от друга? - Да, независимы. Это случайная выборка людей. 2.Выборки независимы друг от друга? - **Нет!**. В группах одни и те же люди (10 человек, каждый пил оба снотворных). **ОСТОРОЖНО**. Эти данные НЕ ПОДХОДЯТ для двухвыборочного t-теста, т.к. группы не являются взаимно независимыми. Нужно использовать **парный** t-тест. <br/> 3.Объем выборки достаточно велик или величины нормально распределены? - Объем выборки мал ```r table(sleep$group) ``` ``` 1 2 10 10 ``` Нужно проверить форму распределения в обеих группах. --- ## Нормально ли распределена зависимая переменная в группах? ```r qqPlot(sleep$extra[sleep$group == 1]) qqPlot(sleep$extra[sleep$group == 2]) ``` <!-- --> Хорошо. Можно аппроксимировать нормальным распределением --- ## Парный t-тест в R Сравним увеличение продолжительности сна при помощи **парного** t-теста. ```r t.test(formula = extra ~ group, data = sleep, paired = TRUE) ``` ``` Paired t-test data: extra by group t = -4.1, df = 9, p-value = 0.003 alternative hypothesis: true mean difference is not equal to 0 95 percent confidence interval: -2.4599 -0.7001 sample estimates: mean difference -1.58 ``` --- ## Опишем результаты ``` Paired t-test data: extra by group t = -4.1, df = 9, p-value = 0.003 alternative hypothesis: true mean difference is not equal to 0 95 percent confidence interval: -2.4599 -0.7001 sample estimates: mean difference -1.58 ``` - Различия изменения продолжительности сна при применении двух препаратов были достоверны ( `\(t_{9} = -4.06\)`, `\(p = <0.01\)`) **Кстати, если бы мы не учли зависимость между группами, мы пришли бы к неверному выводу.** В этом можно убедиться, выполнив этот код: ```r t.test(formula = extra ~ group, data = sleep) ``` --- ## График со средними и доверительными интервалами ```r ggplot(data = sleep, aes(x = factor(group), y = extra)) + stat_summary(fun.data = mean_cl_normal) ``` <!-- --> --- ## Take-home messages - По центральной предельной теореме выборочные средние нормально распределены `$$\bar X \sim N (\mu, \sigma / \sqrt{n})$$`. - Доверительный интервал к среднему значению можно построить из нормального распределения, если известна `\(\sigma\)` в генеральной совокупности, или из `\(t\)`-распределения, если она не известна. - При тестировании нулевой гипотезы оценивают вероятность получения данного или более экстремального значения тестовой статистики при условии, что `\(H_0\)` справедлива. Если эта вероятность меньше выбранного критического уровня значимости, то нулевую гипотезу отвергают. - Для сравнения выборочных средних лучше использовать t-критерий в модификации Уэлча, который учитывает, что дисперсии в группах могут быть разными. - Средние в независимых выборках нужно сравнивать двухвыборочным t-тестом. - Если выборки зависимы — нужно это учесть, поэтому для сравнения средних используется парный t-тест. --- ## Дополнительные ресурсы - OpenIntro: Statistics - Quinn, Keough, 2002 - Sokal, Rohlf, 1995 - Zar, 1999