Задание

• Загрузите объект population. Это вектор, характеризующий размеры особей в некоторой гипотетической популяции.

load("data/population.RData")

• Сделайте из этой совокупности случайную выборку в 5 измерений.

my_sample <-  sample(x = population, size = 5)

my_sample  

[1] 171.6 161.0 166.0 161.7 170.7

• Определите среднее значение и среднеквадратичное отклонение.

mean(my_sample)

[1] 166.2

sd(my_sample)

[1] 4.918

• Проведите такую выборку несколько раз.

Выборочные оценки

• Изучая интересующую нас систему, мы не знаем значений параметров распределения, которое описывает поведение вероятности в генральной совокупности (population).

• Мы проводим оценку (estimate) этих параметров, используя выборку (sample).

• Выборка имеет ограниченный объем(N, или n).

• Выборка попадает в анализ в виде вектора значений \(\{x_1, x_2, x_3, \dots , x_n\}\)

Выборочные оценки

Выборка может быть описана с помощью “5 number statistic”:

Выборочные оценки

Другой способ описания выборки

Как можно судить о свойствах генеральной совокупности по выборке?

Центральная предельная теорема (ЦПТ) говорит, что если мы возьмем достаточно большую выборку из генеральной совокупности, то среднее значение будет нормально распределено с параметрами \(\mu_{\bar x}\) и \(\sigma _{\bar{x}}\):

\[\bar X \sim N (\mu_{\bar x}, \sigma_{\bar x})\]

При чем \(\sigma_{\bar x} = \frac{\sigma}{\sqrt{n}}\).

Важно: это будет так при больших объемах выборки ( \(n > 30\), или даже \(n > 100\)), даже если \(x\) в генеральной совокупности не подчиняется нормальному распределению.

Давайте проверим на опыте, так ли это.

Цена алмазов

Представим, что данные об алмазах из датасета diamonds (пакет gplot2) — это генеральная совокупность.

Перед вами распределение цены алмазов. Давайте будем брать из этого распределения выборки и оценивать по ним среднее значение.

Центральная предельная теорема очень важна в статистике

Если выполняется центральная предельная теорема…

Доверительный интервал
из нормального распределения

Доверительный интервал
из нормального распределения

qnorm(p = c(0.025, 0.975))

[1] -1.96  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-распределение, или распределение Стьюдента

Доверительный интервал из t-распределения

Условия применимости

Выполняются условия, при которых справедлива ЦПТ:

• Наблюдения в выборке независимы друг от друга.
• Большой объем выборки и нет “выбросов” или нормальное распределение \(x\)

Смысл 95% доверительного интервала

Расчет и изображение доверительного интервала в 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

# Границы доверительного интервала
.mean - err

[1] 3826

.mean + err

[1] 4032

Можем записать среднюю цену бриллиантов хорошей огранки и ее доверительный интервал: \(3928.9 \pm 103 (\pm SE \cdot t_{\alpha, df})\)

Строим доверительные интервалы в ggplot

## 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)

Постройте график средней цены с доверительными интервалами для бриллиантов разного качества огранки и прозрачности.

Решение

Вычислим границы доверительного интервала

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

# Границы доверительного интервала
.mean - err

[1] 2086

.mean + err

[1] 2460

Решение: cпособ 1 (некрасивый, требует доработки)

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

ggplot(data = diamonds, aes(x = cut, y = price, colour = clarity)) +
  stat_summary(geom = 'pointrange', fun.data = mean_cl_normal)

Решение: cпособ 1 (доработка)

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

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 (удовлетворительный)

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

На примере доверительных интервалов:

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.

Сравнение выборок

Различия между выборками не всегда видны невооружённым глазом.

tres caracoles by Alberto Villen on Freeimages.com

Нулевая и альтернативная гипотезы

Это первый шаг

tres caracoles by Alberto Villen on Freeimages.com

Нулевая и альтернативная гипотезы — это “два мира”

Вне зависимости от нас, реальность может находиться в одном из двух состояний:

Верные и неверные решения

Тестирование гипотез: Тестовые статистики

Размер кладки черепах

Разберемся с одновыборочным \(t\)-тестом на вымышленном примере.

Представьте, что в одной статье сказано, что средняя плодовитость черепах определенного вида — 8 яиц в кладке.

В вашей выборке из \(35\) черепах — \(\bar x = 9\), \(sd = 2\).

Отличается ли приведенное в статье значение от того, что наблюдается в изученной популяции. Можно ли считать, что авторы статьи изучали что-то иное (другой вид, или популяцию из иных условий).

Gopher Tortoise by Judy Gallagher on Flickr

Одновыборочный t-тест

• \(H_0: \mu = \mu_0\) — Реальная средняя плодовитость черепах такая, как в статье.
• \(H_A: \mu \ne \mu_0\) — Средняя плодовитость отличается от того, что написано в статье.

\(\mu_0\) — это какое-то конкретное значение. В нашей задаче это — 8 яиц в кладке.

\[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\)?

Уровень значимости (p-value)

Уровень значимости (p-value)

2 * ( 1 - pt(2.96, df = 35 - 1) )

[1] 0.005571

Критический уровень значимости

Принимаем решение

Заблуждения о p-values

Гипотезы в двухвыборочном \(t\)-тесте и тестовая статистика

\(H_0: \mu_1 - \mu_2 = 0\) — средние значения не различаются в двух группах

\(H_A: \mu_1 - \mu_2 \ne 0\) — средние значения различаются

Т.е. нас интересует разность выборочных средних.

Ее ожидаемое значение при \(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}}\]

Стандартная ошибка разности средних в 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}}}\]

Приблизительное число степеней свободы

Можно рассчитать по уравнению Уэлча-Саттеруэйта. Это решит проблему.

\[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-теста

Почти такие же, как условия справедливости ЦПТ

• Наблюдения независимы друг от друга.
• Выборки независимы друг от друга (новое условие).
• Объем выборки достаточно велик или величины нормально распределены.

Пример: Гормоны и артериальная гипертензия

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

В датасете Cushings (пакет MASS) записаны данные о секреции двух метаболитов при разных типах синдрома (данные из кн. Aitchison, Dunsmore, 1975).

• Tetrahydrocortisone — секреция тетрагидрокортизона с мочой (мг/сут.)
• Pregnanetriol — секреция прегнантриола с мочой (мг/сут.)
• Type — тип синдрома:
  • a — аденома
  • b — двусторонняя гиперплазия
  • c — карцинома
  • u — не известно

Давайте сравним секрецию тетрагидрокортизона при аденоме и двусторонней гиперплазии надпочечников. Различается ли она?

Открываем данные

library(MASS)
data("Cushings")

Все ли правильно открылось?

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

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 ...

Знакомимся с данными

Есть ли пропущенные значения?

colSums(is.na(Cushings))

Tetrahydrocortisone       Pregnanetriol                Type 
                  0                   0                   0 

Каковы объемы выборки в каждой группе?

table(Cushings$Type)

 a  b  c  u 
 6 10  5  6 

Обратите внимание, объемы выборок маленькие.

Проверяем условия применимости…

1.Наблюдения независимы друг от друга?

• Да, независимы. Это случайная выборка.

2.Выборки независимы друг от друга?

• Да, независимы. В группах разные люди (естественно, т.к. тип синдрома у человека может быть только какой-то один).

3.Объем выборки достаточно велик или величины нормально распределены?

• Объем выборки мал

table(Cushings$Type)

 a  b  c  u 
 6 10  5  6 

Нужно проверить форму распределения в обеих группах.

Нормально ли распределены концентрации тетрагидрокортизона в группах?

library(car)
qqPlot(Cushings$Tetrahydrocortisone[Cushings$Type == 'a'])
qqPlot(Cushings$Tetrahydrocortisone[Cushings$Type == 'b'])

Удовлетворительно. При таких малых объемах выборки сложно ожидать лучшего.

Двухвыборочный t-тест в R: способ 1.

Сравним секрецию тетрагидрокортизона при помощи двухвыборочного t-теста.

t.test(x = значения_в_гр.1, 
       y = значения_в_гр.2)

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\))

Можно указать в скобках не сравнение с \(\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 = логический_вектор)

t.test(formula = Tetrahydrocortisone ~ Type, data = Cushings, 
       subset = Cushings$Type %in% c('a', 'b'))

Результаты точно такие же, естественно.

Задание 4

Посмотрите структуру результатов (tt) при помощи функции str() и извлеките из них:

• степени свободы
• уровень значимости
• значение t-критерия

Решение: Что спрятано в результатах?

Как называются отдельные элементы результатов можно узнать посмотрев их структуру при помощи функции str()

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"

Решение: Можно получить элементы результатов в виде отдельных чисел

tt$parameter # степени свободы

   df 
10.69 

tt$p.value # уровень значимости

[1] 0.001719

tt$statistic # значение t-критерия

    t 
-4.15 

График со средними и доверительными интервалами

ggplot(data = Cushings[Cushings$Type %in% c('a', 'b'), ],
       aes(x = Type, y = Tetrahydrocortisone)) +
  stat_summary(fun.data = mean_cl_normal)

Подготовка данных

a <- Cushings$Tetrahydrocortisone[Cushings$Type == 'a']
b <- Cushings$Tetrahydrocortisone[Cushings$Type == 'b']

mean_a <- mean(a)
sd_a <- sd(a)
n_a <- length(a)
se_a <- sd_a/sqrt(n_a)

mean_b <- mean(b)
sd_b <- sd(b)
n_b <- length(b)
se_b <- sd_b/sqrt(n_b)

t_emp <- abs(mean_a - mean_b)/sqrt(se_a^2 + se_b^2)

Пермутация векторов

all_data <- sample(c(a, b))
a_perm <- all_data[1:n_a]
b_perm <- all_data[(n_a + 1):length(all_data)]

mean_a <- mean(a_perm)
sd_a <- sd(a_perm)
n_a <- length(a_perm)
se_a <- sd_a/sqrt(n_a)

mean_b <- mean(b_perm)
sd_b <- sd(b_perm)
n_b <- length(b_perm)
se_b <- sd_b/sqrt(n_b)

t_perm <- abs(mean_a - mean_b)/sqrt(se_a^2 + se_b^2)

Многократные пермутации

t_perm_vector <- rep(NA, 1000)

for(i in 1:999){
    all_data <- sample(c(a, b))
    a_perm <- all_data[1:n_a]
    b_perm <- all_data[(n_a +1):length(all_data)]

    mean_a <- mean(a_perm)
    sd_a <- sd(a_perm)
    n_a <- length(a_perm)
    se_a <- sd_a/sqrt(n_a)

    mean_b <- mean(b_perm)
    sd_b <- sd(b_perm)
    n_b <- length(b_perm)
    se_b <- sd_b/sqrt(n_b)

    t_perm <- abs(mean_a - mean_b)/sqrt(se_a^2 + se_b^2)

    t_perm_vector[i] <-  t_perm
}


  t_perm_vector[1000] <- t_emp

Отображаем графически

t_perm_vector_dataframe <- as.data.frame(t_perm_vector)
ggplot(data=t_perm_vector_dataframe, aes(x=t_perm_vector))+geom_histogram()+geom_vline(xintercept = t_emp, colour = 'red', size = 2)

Как часто мы встречаем значение t-статистики как для изначальных выборок или больше?

Результаты параметрического теста

mean(t_perm_vector >= t_emp) #Доля пермутационных статистик которые больше или равны выборочной статистике

[1] 0.001

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 

Задание 5

Файл aml.csv содержит данные о влиянии регулярной химиотерапии на продолжительность ремиссии.

Прочитаем эти данные

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-теста.

Решение: Разведочный анализ

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 ...

# Делаем group фактором
rem$group <- factor(rem$group, labels = c("therapy", "control"))

# Пропущенные значения?
colSums(is.na(rem))

 time group 
    0     0 

# Объемы выборок?
table(rem$group)

therapy control 
     11      12 

Решение: Проверка условий применимости

# Выборки малы, поэтому проверяем на нормальность и отсутствие выбросов
op <- par(mfrow = c(1, 2))
qqPlot(rem$time[rem$group == "therapy"])

[1] 11  1

qqPlot(rem$time[rem$group == "control"])

[1] 12 11

par(op)

Решение: Проверка условий применимости

Есть отскакивающее значение в группе therapy, в этой группе оно 11 по порядку (видно на qqPlot()).

rem$time[rem$group == "therapy"]

 [1]   9  13  13  18  23  28  31  34  45  48 161

rem$time[rem$group == "therapy"][11]

[1] 161

Пусть в этом случае это вполне реальное наблюдение, но это значение придется исключить, т.к. t-тест не устойчив к выбросам.

Решение: Двухвыборочный t-тест

# удалено наблюдение со значением 161
tt <- t.test(formula = time ~ group, data = rem,
       subset = rem$time != 161)

Или то же самое:

# удалено 11 наблюдение в группе `therapy`
tt <- t.test(x = rem$time[rem$group == 'therapy'][-11],
       y = rem$time[rem$group == 'control'])

Решение: Результаты двухвыборочного t-теста

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\)).

Решение: График со средними и доверительными интервалами

ggplot(data = rem[rem$time != 161, ], aes(x = factor(group), y = time)) + 
  stat_summary(fun.data = mean_cl_normal)

Пример: Снотворное

В датасете sleep содержатся данные об увеличении продолжительности сна по сравнению с контролем после применения двух снотворных препаратов (Cushny, Peebles, 1905, Student, 1908)

Одинаково ли два снотворных влияют на увеличение продолжительности сна?

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

# str(sleep)

Проверяем условия применимости…

1. Наблюдения независимы друг от друга?

• Да, независимы. Это случайная выборка людей.

2. Выборки независимы друг от друга?

• Нет!. В группах одни и те же люди (10 человек, каждый пил оба снотворных).
  ОСТОРОЖНО. Эти данные НЕ ПОДХОДЯТ для двухвыборочного t-теста, т.к. группы не являются взаимно независимыми. Нужно использовать парный t-тест.

3. Объем выборки достаточно велик или величины нормально распределены?

• Объем выборки мал

table(sleep$group)

 1  2 
10 10 

Нужно проверить форму распределения в обеих группах.

Нормально ли распределена зависимая переменная в группах?

qqPlot(sleep$extra[sleep$group == 1])
qqPlot(sleep$extra[sleep$group == 2])

Хорошо. Можно аппроксимировать нормальным распределением

Парный t-тест в R

Сравним увеличение продолжительности сна при помощи парного t-теста.

t.test(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2], paired = TRUE)

    Paired t-test

data:  sleep$extra[sleep$group == 1] and sleep$extra[sleep$group == 2]
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:  sleep$extra[sleep$group == 1] and sleep$extra[sleep$group == 2]
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\))

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

t.test(formula = extra ~ group, data = sleep)

График со средними и доверительными интервалами

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