Статистика и R
Занятие 8
10 апреля 2023
Занятие 8
Статистика
- ANOVA
- Непараметрические аналоги ANOVA
- Проблема множественного тестирования
R
- R Markdown
- Quarto
R Markdown
Rmd документ
File -> New File -> R Markdown ...
Части Rmd документа
- YAML шапка [YAML header] - отделяется
---и содержит информацию о типе генерируемого документа - Блок кода [code chunk] - содержит код
- Заголовки - выделяются одним или несколькими символами
# - Обычный текст
Создать отчет из Rmd документа
- Можно работать как с R скриптом либо
Knit -> Knit to HTML
Markdown разметка
Заголовки
# Заголовок первого порядка
## Заголовок второго порядка
### Заголовок третьего порядка
Списки
- Такой элемент списка
- Другой элемент списка
- Вот еще такой элемент
- И такой
- Первый элемент списка
- Второй элемент списка
+ Пример 1
+ Пример 2
Опции блока кода
eval=FALSE- не выполнять код (но показать его в отчете)echo=FALSE- не показывать код (но показать результат его выполнения)include=FALSE- выполнить код, но не включать в отчет ни код, ни его реузльтатerror=TRUE- допускать ошибки в коде
Сравнение трех и более выборок
ANOVA
- ANalysis Of VAriance
- тест с одной гипотезой
- тест на равенство средних среди нескольких групп
Варианты ANOVA
Однофакторный дисперсионный анализ [One-Way ANOVA] - один фактор (категориальная независимая переменная) и несколько уровней фактора.
Двухфакторный дисперсионный анализ [Two-Way ANOVA] - два фактора и несколько уровней факторов.
Многофакторный дисперсионный анализ [Factorial ANOVA]
Требования ANOVA
- случайные независимые выборки (наблюдения независимы)
- данные в выборках распределены приблизительно нормально
- дисперсия данных в выборках одинаковая (гомоскедастичность)
Нормальность
Гомоскедастичность
- Оценить по выборочному стандартному отклонению
Гомоскедастичность
- Посмотреть на график остатков (residuals) - xi−¯x
Гомоскедастичность
- Тест Левена на равенство дисперсий (H0: дисперсии равны)
А если дисперсии не равны?
Если выборочные дисперсии сильно отличаются (часто если размеры выборок сильно отличаются), то надежнее использовать Whelch One-Way ANOVA - этот тест не трубует гомоскедастичности.
Однофакторный ANOVA
- H0:μ1=μ2=...=μk - все средние равны
- Хотя бы одно среднее отличается
При чем здесь дисперсия?
F=S2betweenS2within, где
- S2within - дисперсия внутри выборок
- S2between - дисперсия между средними выборок
Если F < 1, это значит, что большой разницы между средими выборок нет.
F-распределение
Чем больше отношение дисперсии между средними выборок к дисперсии внутри выборок, тем больше значение F-статистики. Большие значения F-статистики говорят о несостоятельности нулевой гипотезы. Для подсчета p-значения используется правый (верхний) хвост F-распределения.
ANOVA - шаг 1
1. Дисперсия между средними выборок (MSG - mean square between groups)
S2between=MSG=1k−1∑ki=1ni(¯xi−¯x)2=1dfGSSG
k - количество сравниваемых выборок/групп,
ni - размер выборки i,
¯xi - среднее выборки i,
¯x - среднее по всем выборкам/группам
SSG - sum of squares between groups,
dfG=k−1
ANOVA - шаг 2
2. Дисперсия внутри выборок (MSE - mean square error)
S2within=MSE=1n−k∑ki=1(ni−1)s2i=1dfESSE
n - общее число наблюдений,
s2i - дисперсия выборки i
SSE - sum of squared errors,
dfE=n−k
ANOVA - шаг 3
3. F-статистика
F=MSGMSE∼F(dfG,dfE)
Пример 1
Два разных удобрения были испытаны на полях с растениями. Есть ли различия в объемах собранной растительной массы на полях, где использовались удобрения, и контрольном поле?
# A tibble: 10 × 4
id ctrl trt1 trt2
<int> <dbl> <dbl> <dbl>
1 1 4.17 4.81 6.31
2 2 5.58 4.17 5.12
3 3 5.18 4.41 5.54
4 4 6.11 3.59 5.5
5 5 4.5 5.87 5.37
6 6 4.61 3.83 5.29
7 7 5.17 6.03 4.92
8 8 4.53 4.89 6.15
9 9 5.33 4.32 5.8
10 10 5.14 4.69 5.26Пример 1
Пример 1
1.
MSG=1k−1∑ki=1ni(¯xi−¯x)2 MSG=13−1(10×(5.03−5.07)2+10×(4.66−5.07)2+ 10×(5.53−5.07)2)
MSG=12×3.813=1.9065
2.
MSE=1n−k∑ki=1(ni−1)s2i MSE=130−3((10−1)×0.34+(10−1)×0.63+(10−1)×0.196)
MSE=127×10.494=0.3887
Пример 1
3.
F=MSGMSE=1.90650.3887=4.9048
p−value= 0.015
Если p−значение<α, то имеет смысл смотреть на попарные сравнения выброк с помощью апостериорных тестов (post-hoc tests).
ANOVA в R
- формула:
количественная переменная ~ качественная переменная - данные: датафрейм
Whelch One-Way ANOVA
One-way analysis of means (not assuming equal variances)
data: weight and group
F = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.01739“Обычный” ANOVA:
H-критерий Краскела-Уоллиса
- непараметрический аналог однофакторного ANOVA - Kruskal-Wallis test
- или расширение критерия Манна-Уитни для более чем двух групп
H-критерий Краскела-Уоллиса
Алгоритм:
- X1,...,Xn,Y1,...,Ym,...,Z1,...,Zl - k независимых выборок
- приписываем ранги всем наблюдениям из объединенных выборок
- считаем суммы рангов наблюдений для каждой выборки: Ri
- H=12N(N+1)∑R2ini−3(N+1), где N - суммарный объем выборок
- H∼χ2(k−1), где k - количество групп
Пример 2
Вы сравниваете 3 географически изолированные популяции амёб по количеству ложноножек. Можно ли на основании собранных данных утверждать (на уровне значимости 0.05), что этих амёб можно разделить хотя бы на два вида, различающихся по количеству ложноножек?
Пример 2
- R1=7+7+9+10+12+14.5=59.5
- R2=2.5+4.5+4.5+12=23.5
- R3=1+2.5+7+12+14.5=37
H=12N(N+1)∑R2ini−3(N+1)
H=1215×(15+1)(59.526+23.524+3725)−3(15+1)
H=0.05×(590.05+138.06+273.8)−48=2.09
df=k−1=2
χ2(2)0.9<2.09<χ2(2)0.1, то есть p-значение лежит в интервале от 0.1 до 0.9. На уровне значимости 0.05 H0 отвергнуть нельзя.
Парные выборки
Параметрические тесты:
- ANOVA для повторных измерений (Repeated measures analysis of variance (rANOVA))
Непараметрические тесты:
- тест Фридмана (Friedman test)
ANOVA на графике
Тест Краскела-Уоллиса на графике
Апостериорные тесты
Зачем нужны post-hoc тесты?
понять, какие именно выборки отличаются
- Можно использовать попарные t-тесты или тесты Манна-Уитни?
При проведении нескольких тестов возникает проблема множественного тестирования.
Одна гипотеза, разные ошибки
Нулевая гипотеза H0 - наши априорные представления о картине мира. Увидев некоторые данные, мы принимаем решение принять или отвергнуть H0.
- ошибка I рода (False Positive) α - H0 верна, но мы ее отвергли,
- ошибка II рода (False Negative) β - H0 неверна, но мы ее приняли.
Мы хотим:
- гарантировать, что вероятность ошибки I рода ниже заранее установленного уровня значимости α (обычно 5%),
- максимизировать мощность критерия 1−β, то есть минимизировать вероятность ошибки II рода.
Ошибки I и II рода
Несколько гипотез, еще больше ошибок
- ошибка I рода для индивидуальной гипотезы
- вероятность хотя бы одной ошибки I рода в m тестах: 1−(1−α)m - это Family-wise Error Rate [FWER]
- частота ошибок I рода (False Positive) среди всех отвергнутых H0: FDR=FPFP+TP - False Discovery Rate [FDR]
| Число групп | Число тестов | FWER |
|---|---|---|
| 2 | 1 | 0.050 |
| 3 | 3 | 0.143 |
| 4 | 6 | 0.265 |
| 5 | 10 | 0.401 |
| 6 | 15 | 0.537 |
| 7 | 21 | 0.659 |
| 8 | 28 | 0.762 |
| 9 | 36 | 0.842 |
| 10 | 45 | 0.901 |
Поправка Бонферрони
Контролирует FWER на уровне α.
Умножаем p-значение (то же, что и делим уровень значимости) на число тестов. Если после такой процедуры p-значение останется достаточно маленьким, то мы продолжаем отвергать нулевую гипотезу.
Допустим мы провели 10000 тестов, в одном из них получили p-значение 2×10−5. Это достаточно низкое p-значение, чтобы отвергнуть H0 при уровне значимости α=0.05?
Padj=2×10−5×104=0.2>α
Поправка Холма-Бонферрони
- сортируем p-значения так, чтобы p1≤p2≤...≤pn
- находим максимальное k такое, что pk≤αn−k+1
- для всех p1,...,pk отвергаем H0
Поправка Шидака-Холма
- сортируем p-значения так, чтобы p1≤p2≤...≤pn
- находим максимальное k такое, что pk≤1−(1−α)1n−k+1
- для всех p1,...,pk отвергаем H0
Поправка Бенджамини-Хохберга
Контролирует FDR на уровне α.
- сортируем p-значения так, чтобы p1≤p2≤...≤pn
- находим максимальное k такое, что pk≤k×αn
- для всех p1,...,pk отвергаем H0
Это более мягкая поправка по сравнению с поправкой Бонферрони, ее стоит использовать, когда действительно много тестов.
Пример 3
Допустим, вы провели 10 тестов и получили следующие p-значения: 0.0001, 0.007, 0.01, 0.05, 0.06, 0.12, 0.25, 0.37, 0.62, 0.83.
Для скольки тестов мы можем отвергнуть нулевые гипотезы (на уровне значимости 0.05) после применения поправок Бонферрони и Бенджамини-Хохберга?
Пример 3
Бонферрони
- нужно проверить pi×n<α
- поправленные p-значения: 0.001, 0.07, 0.1, 0.5, 0.6, 1.0, 1.0, 1.0, 1.0, 1.0
Бенджамини-Хохберг
- αn=0.005
- для k=1, p1<0.005
- для k=2, p2<0.01
- для k=3, p3<0.015
- для k=4, p4>0.02
Пример 3
Бонферрони
pvals <- c(0.0001, 0.007, 0.01, 0.05, 0.06, 0.12, 0.25, 0.37, 0.62, 0.83)
p.adjust(pvals, method = "bonferroni") [1] 0.001 0.070 0.100 0.500 0.600 1.000 1.000 1.000 1.000 1.000Бенджамини-Хохберг
Тест Тьюки
Tukey Honest Significant Differences [Tukey HSD]
- после ANOVA (использует внутреннюю поправку на множественное тестирование, контролирующую FWER)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = weight ~ group, data = PlantGrowth)
$group
diff lwr upr p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1 0.865 0.1737839 1.5562161 0.0120064diff- разница между средними,lwr, upr- нижняя и верхняя граница 95% (по умолчанию) доверительного интервала,p adj- p-value после поправки на множественное тестирование.
Тест Данна
Dunn test
- после теста Краскела-Уоллиса
Kruskal-Wallis rank sum test
data: weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842# A tibble: 3 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 weight ctrl trt1 10 10 -1.12 0.264 0.264 ns
2 weight ctrl trt2 10 10 1.69 0.0912 0.137 ns
3 weight trt1 trt2 10 10 2.81 0.00500 0.0150 * Тест Даннетта
Dunnett test
- после ANOVA
- интересуют только сравнения некоторых групп с “контрольной” группой.
# не делает сам поправку на МТ
DescTools::DunnettTest(weight ~ group, data = PlantGrowth, control = "ctrl")
Dunnett's test for comparing several treatments with a control :
95% family-wise confidence level
$ctrl
diff lwr.ci upr.ci pval
trt1-ctrl -0.371 -1.0215695 0.2795695 0.3227
trt2-ctrl 0.494 -0.1565695 1.1445695 0.1535
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Тест Даннетта
Поправка на множественное тестирование.
Тест Геймса-Ховелла
Games-Howell test
- после Welch ANOVA (если дисперсии не равны).
# A tibble: 3 × 8
.y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
* <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 weight ctrl trt1 -0.371 -1.17 0.430 0.475 ns
2 weight ctrl trt2 0.494 -0.101 1.09 0.113 ns
3 weight trt1 trt2 0.865 0.114 1.62 0.024 *