Первый шаг в сравнении – формулировка нулевой гипотезы.

Положение H0 привилегированное: именно ее мы тестируем, принимаем или отвергаем.

Чаще всего она формулируется как отсутствие различий между сравниваемыми объектами. Например: ящерицы из двух популяций одинаково крупные.

Вместе с нулевой гипотезой рождается на свет и альтернативная гипотеза.

В общем виде, она формулируется как присутствие различий и включает все частные случаи, например “размеры ящериц из двух популяций неодинаковы”.

Решение о том, принять или отвергнуть H0 мы принимаем после статистического теста. Но результат теста не тождественен реальному состоянию дел - он не скажет, верна ли H0, только укажет, как с ней поступить.

Вне зависимости от нас, реальность может находиться в одном из двух состояний: H0 верна, ящерицы одинаковы, либо H0 неверна, и ящерицы различаются.

Таким образом, возможно четыре исхода теста.

В мире где ящерицы одинаковы мы можем

• принять H0 (и это будет верное решение),
• либо отвергнуть ее (и это будет ошибка).

Аналогично, в мире где ящерицы различаются, мы можем

• принять H0 (что ошибочно),
• либо отвергнуть ее (что верно).

Ошибка I рода: нашли то, чего нет

Ошибка II рода: пропустили то, что было

Вероятность ошибки I рода мы задаем сами — это уровень значимости теста, \(\alpha\). Это та вероятность, меньше которой мы отказываемся верить в справедливость H0.

Ошибка второго рода обычно скрыта. Но ее можно уменьшить вслепую, например увеличив объем выборки.

Статистические таблицы

Принцип прост – мы вычисляем эмпирическое значение тестовой статистики, и сравниваем его с критическим значением этой же статистики из таблицы.

Т.е. мы сравниваем тестовую статистику с тестовой статистикой: t-критерий с t-критерием, \(\chi^2\) с \(\chi^2\).

P-value

Принятие решения на основании p-value стало возможным когда стали доступны большие вычислительные мощности.

P-value отражает вероятность обнаружить тот уровень различий между выборками который мы видим, при условии что H0 верна.

Если р велико, мы решаем что выборки принципиально одинаковы, если мало - что выборки принципиально различаются.

Принимая решение на основании p-value мы сравниваем ее с уровнем значимости \(\alpha\), т.е. сравниваем вероятность с вероятностью.

Допустим, мы сравнили выборки и получили р = 0.03

P-value = 0.03 не значит, что H0 верна с вероятностью 3%!

P-value = 0.03 значит, что в мире, где выборки одинаковы а H0 верна, шанс получить результат который мы получили составляет 3%.

Уже мы сами решаем, кажется ли нам такая вероятность приемлемой.

Из эмпирического распределения X создадим распределение Z, где каждое значение \(x_i\) будет заменено на \(z_i\).

1. Выполним центрирование. Для этого заменим значения переменной на девиаты: \(x_i - \bar{x}\)

2. Нормируем разброс. Для этого разделим девиаты на стандартное отклонение выборки.

\[z_i=\frac{x_i - \bar{x}}{SD}\]

Давайте познакомимся с тремя полезными функциями для работы с распределениями: r, q, p

r = random number generation.

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

Например: rnorm(1000, mean = 20, sd = 2) сгенерирует выборку в 1000 случайных значений из нормального распределения.

q = quantile function.

С ее помощью можно получить квантиль (точнее персентиль) - то значение переменной, которое отсекает заданную часть распределения.

Например: qnorm(0.5, mean = 20, sd = 2) вернет среднюю, т.к. ровно 50% значений в нормальном распределении \(< \mu\)

qnorm(0.025, mean = 20, sd = 2) вернет то значение X, которое делит распределение на куски в 2.5% и 97.5%, т.е. 2.5% всех значений будут меньше, а 97.5% - больше него.

С помощью функции qnorm получите 5-number summary для распределения с параметрами \(\mu\) = 3 и \(\sigma\) = 5

Для 5-number summary нам нужны минимальное и максимальное значение, медиана, I и III квантили.

Создадим вектор

five_numbers <- c(0, 0.25, 0.5, 0.75, 1)

и передадим его функции qnorm

qnorm(five_numbers, 3, 5)

## [1]   -Inf -0.372  3.000  6.372    Inf

p = probability distribution function

С ее помощью можно рассчитать вероятность того, что случайная величина, взятая из данного распределения, окажется меньше заданного нами значения.

Эта операция принципиально обратна тому, что делает q.

Например: pnorm(20, mean = 20, sd = 2) вернет вероятность 0.5, или 50%, поскольку мы передали ей в качестве аргумента среднюю

Эта функция работает кумулятивно:

round(pnorm(c(0, 15, 20, 25, 40), 20, 2), 3)

## [1] 0.000 0.006 0.500 0.994 1.000

\[t=\frac{d}{SE_d}\]

• \(d=\bar{x_1} - \bar{x_2}\) - это разность между двумя средними значениями

• \(SE_d\) - Общее среднеквадратичное отклонение разности двух средних

\[SE_d = \sqrt{\frac{sd_1^2(n_1-1) +sd_2^2(n_2-1)}{n_1+n_2-2}(\frac{1}{n_1} + \frac{1}{n_2})}\]

Если \(n_1 = n_2\), то формула существенно упрощается

\[SE_d = \sqrt{\frac {sd_1^2} {n_1} + \frac {sd_2^2} {n_2}}\]

Таким образом, t-распределение это стандартизованное распределение разностей двух средних значений из одной генеральной совокупности

Распределение t-статистики описывает заковыристая функция. Нам нужно знать про нее один важный факт: форма распределения зависит от единственного параметра \(df\) - числа степеней свободы.

\[df = n_1 + n_2 - 2\]

Давайте посмотрим, как выглядят распределения с разными df:

Давайте выполним t-тест и решим, как нам поступить с нулевой гипотезой (для этого нам пригодятся операции qt и pt)

Для начала создадим две выборки длин ящериц из популяций Берлина и Саратова. В этих гипотетических выборках длины распределены нормально, и имеют заведомо отличающимися \(\mu\)

# Зерно для генератора случайных чисел для сопоставимости результатов
set.seed(456) 
# Создаем две выборки по 100 из нормального распределения с разными параметрами
Saratov <- rnorm(n = 100, mean = 130, sd = 5)
Berlin <- rnorm(n = 100, mean = 129, sd = 5)
city <- c(rep("B", 100), rep("S", 100))
# Сохраняем выборки в датафрейме для удобства
lizards <- data.frame(city = factor(city),
                        length = c(Berlin, Saratov))


head(lizards)

##   city length
## 1    B    130
## 2    B    133
## 3    B    129
## 4    B    129
## 5    B    121
## 6    B    135

library(ggplot2) # Загрузим библиотеку
theme_set(theme_bw()) # Зададим тему

Сконструируем “скелет” графика

ggplot(lizards, aes(x = length)) + 
  geom_histogram()

Здесь ящерицы из разных мест пока еще смешаны.

ggplot(lizards, aes(x = length)) + 
  geom_histogram(binwidth = 3)

gg_length <- ggplot(lizards, aes(x = length, fill = city)) + 
    geom_histogram(binwidth = 3, colour = "grey40", position = "dodge")
gg_length

gg_length + 
  labs(x = "Length (cm)", 
       y = "Count", 
       title ="Length distribution of lizards", 
       fill = "City")

Наш график готов!

t_lizards <- t.test(data = lizards, length ~ city)
# length ~ city показывает, что значения переменной
# length сгруппированы по признаку city
# Порядок записи важен!
t_lizards

## 
##  Welch Two Sample t-test
## 
## data:  length by city
## t = -3, df = 200, p-value = 0.002
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.562 -0.811
## sample estimates:
## mean in group B mean in group S 
##             128             131

t_lizards это комплексный объект, который можно изучить с помощью функции str()

Нас интересуют две переменные: значение тестовой t-статистики

t_value_lizards <- t_lizards$statistic
t_value_lizards

##     t 
## -3.13

и итоговое p-value

p_value_lizards <- t_lizards$p.value
p_value_lizards

## [1] 0.00199

Рассчитаем “табличное” значение t-критерия с помощью функции qt()

• так как критерий двухсторонний, мы должны разбить 5% на два кусочка по 2.5%, т.е. p = c(0.025, 0.975)
• число степеней свободы \(n_1 + n_2 - 2\) равно 100 + 100 - 2 = 198

Подставим эти аргументы:

qt(p = c(0.025, 0.975), 198)

## [1] -1.97  1.97

Теперь с этими критическими значениями можно сравнить эмпирический результат:

t_value_lizards

##     t 
## -3.13

С помощью команды pt() определим вероятность того, что случайная величина, взятая из t-распределения, окажется меньше рассчитанной нами t-статистики

pt(t_value_lizards, df = 198)

##        t 
## 0.000994

Это “сырая” величина. Её еще нельзя сравнивать с уровнем значимости.

Мы должны либо умножить p.value на 2, либо поделить \(\alpha\) пополам, и только после этого сравнивать их.

p.value <- 2 * pt(t_value_lizards, df = 198)
p.value

##       t 
## 0.00199

Этот результат совпадает с тем, что рассчитал t-тест

Вопрос: Вероятность какого события отражает p=0.002?

Это вероятность того, что выборки которые мы имеем были получены из одной совокупности.

Иными словами:

Это НЕ вероятность того, что H0 верна.

Это вероятность того, что в мире, где H0 верна, а длины ящериц из Берлина и Саратова равны, шанс получить эмпирические выборки с такими различиями как у нас, составляет 0.2%.

Кстати, какова разница между средними длинами в наших выборках?

• Независимость выборок друг от друга
  
• Независимость наблюдений внутри групп
• Нормальное распределение сравниваемых величин (если выборки малы, n < 30)
  
• Равенство дисперсий (можно нарушать, требуется коррекция по методу Велча)

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

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

rem <- read.csv("data/2.2_aml.csv", header = 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 ...

• В переменной time представлена продолжительность ремиссии в днях.
• group указывает, к какой экспериментальной группе принадлежал пациент. В группе 1 проводилась регулярная химиотерапия, в группе 2 - нет.

Ваша задача сравнить эти группы с помощью t-теста.

library(car)
par(mfrow = c(1, 2))
qqPlot(rem$time[rem$group == 1])

## [1] 11  1

qqPlot(rem$time[rem$group == 2])

## [1] 12 11

par(mfrow = c(1, 1))

Есть выброс, хорошо бы проверить это наблюдение.

rem1 <- rem[rem$time < 150, ] # исключаем выброс
t.test(data = rem1, time ~ group)

## 
##  Welch Two Sample t-test
## 
## data:  time by group
## t = 0.8, df = 20, p-value = 0.4
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.49 17.39
## sample estimates:
## mean in group 1 mean in group 2 
##            26.2            21.2

Или так:

t.test(rem1$time[rem1$group==1], rem1$time[rem1$group==2])

Статистическая значимость бывает только одна и она говорит о том, в каких отношениях состоят обнаруженные нами эмпирически различия с \(\alpha\).

Если мы отвергаем H0, имеются различия.

Если мы не отвергаем H0, все одинаково.

Пожалуйста, никаких недостоверных различий! Если различия есть, то они достоверны. Если различия не достоверны, то их нет.

А что делать, если сердцем чувствуешь, что что-то должно быть, но тест не дает нужного ответа? Собрать выборки побольше, и посчитать заново. С уже проведенным тестом спорить не надо.