class: middle, left, inverse, title-slide .title[ # Линейные модели с дискретными предикторами ] .subtitle[ ## Линейные модели… ] .author[ ### Марина Варфоломеева, Вадим Хайтов ] .date[ ### Осень 2022 ] --- # Линейные модели с дискретными предикторами (дисперсионный анализ) ## Вы сможете - Объяснить, в чем опасность множественных сравнений, и как с ними можно бороться - Рассказать, как в дисперсионном анализе моделируются значения зависимой переменной - Интерпретировать и описать результаты, записанные в таблице дисперсионного анализа - Перечислить и проверить условия применимости дисперсионного анализа - Провести множественные попарные сравнения при помощи post hoc теста Тьюки, представить и описать их результаты - Построить график результатов дисперсионного анализа --- ## Дисперсионный анализ (Analysis Of Variance, ANOVA) __Дисперсионный анализ в широком смысле__ --- анализ изменений непрерывной зависимой переменной в связи с разными источниками изменчивости (предикторами). Мы использовали его для тестирования значимости предикторов в линейных моделях. -- __Дисперсионный анализ в узком смысле__ --- это частный случай, когда в линейной модели используются только дискретные предикторы (факторы). Он используется для сравнения средних значений зависимой переменной в дискретных группах, заданных факторами.. --- ## Пример: яйца кукушек Различаются ли размеры яиц кукушек в гнездах разных птиц-хозяев? Датасет `cuckoos` из пакета `DAAG`: - `species` --- вид птиц-хозяев (фактор) - `length` --- длина яиц кукушек в гнездах хозяев (зависимая переменная) .small[Данные: Latter, 1902; источник: Tippett, 1931] --- ## Открываем данные ```r library(DAAG) data("cuckoos") # Положим данные в переменную с коротким названием, чтобы меньше печатать eggs <- cuckoos head(eggs, 3) ``` ``` length breadth species id 1 21.7 16.1 meadow.pipit 21 2 22.6 17.0 meadow.pipit 22 3 20.9 16.2 meadow.pipit 23 ``` ```r # Сократим названия переменных colnames(eggs) <- c('len', 'br', 'sp', 'id') ``` --- ## Изменим названия уровней фактора, чтобы было легче понять о каких птицах речь ```r levels(eggs$sp) ``` ``` [1] "hedge.sparrow" "meadow.pipit" "pied.wagtail" "robin" "tree.pipit" [6] "wren" ``` ```r levels(eggs$sp) <- c("ЛесЗав", "ЛугКон", "БелТряс", "Малин", "ЛесКон", "Крапив") ``` --- ## Исследуем данные ```r # Пропущенных значений нет colSums(is.na(eggs)) ``` ``` len br sp id 0 0 0 0 ``` ```r # Данные не сбалансированы (размеры групп разные) table(eggs$sp) ``` ``` ЛесЗав ЛугКон БелТряс Малин ЛесКон Крапив 14 45 15 16 15 15 ``` --- ## Задание 1 Дополните код, чтобы построить график зависимости размера яиц кукушек (`len`) от вида птиц-хозяев (`sp`), в гнездах которых были обнаружены яйца. На графике должны быть изображены средние значения и их 95% доверительные интервалы, а цвет должен соответствовать виду птиц-хозяев. ```r library() theme_set( ) # График средних с 95% доверительными интервалами ggplot(data = , aes()) + stat_summary(geom = , fun.data = ) ``` <!-- --> --- ## Решение 1 ```r library(ggplot2) theme_set(theme_bw()) ``` ```r # График средних с 95% доверительными интервалами ggplot(data = eggs, aes(x = sp, y = len, colour = sp)) + stat_summary(geom = "pointrange", fun.data = mean_cl_normal) ``` <!-- --> --- ## "Некрасивый" порядок уровней на графике На этом графике некрасивый порядок уровней: средние для разных уровней фактора `eggs$sp` расположены, как кажется, хаотично. Порядок групп на графике определяется порядком уровней фактора. ```r # "старый" порядок уровней levels(eggs$sp) ``` ``` [1] "ЛесЗав" "ЛугКон" "БелТряс" "Малин" "ЛесКон" "Крапив" ``` --- ## Меняем порядок уровней Давайте изменим порядок уровней в факторе `eggs$sp` так, чтобы он соответствовал возрастанию средних значений длины яиц `eggs$len`. ```r # "старый" порядок уровней levels(eggs$sp) ``` ``` [1] "ЛесЗав" "ЛугКон" "БелТряс" "Малин" "ЛесКон" "Крапив" ``` ```r # переставляем уровни в порядке следования средних значений eggs$sp <- reorder(eggs$sp, eggs$len, FUN = mean) # "новый" порядок уровней стал таким levels(eggs$sp) ``` ``` [1] "Крапив" "ЛугКон" "Малин" "БелТряс" "ЛесКон" "ЛесЗав" ``` --- ## График с новым порядком уровней С новым порядком уровней нам легче визуально сравнивать друг с другом категории. Поскольку, изменив порядок уровней, мы внесли изменения в исходные данные, придется полностью обновить график (т.к.`ggplot()` хранит данные внутри графика). ```r ggplot(data = eggs, aes(x = sp, y = len, colour = sp)) + stat_summary(geom = "pointrange", fun.data = mean_cl_normal) ``` <!-- --> --- ## Понравившийся график, если понадобится, можно в любой момент довести до ума, а остальные удалить ```r gg_means <- ggplot(data = eggs, aes(x = sp, y = len, colour = sp)) + stat_summary(geom = "pointrange", fun.data = mean_cl_normal, size = 1) + labs(x = "Вид хозяев", y = "Длина яиц кукушек, мм") + scale_colour_brewer(name = "Вид \nхозяев", palette = "Dark2") + theme(legend.position = "none") gg_means ``` <!-- --> --- class: middle, center, inverse # Множественные сравнения --- ## Множественные сравнения Мы могли бы сравнить длину яиц в гнездах резных хозяев при помощи t-критерия. У нас всего 6 групп. Сколько возможно между ними попарных сравнений? .pull-left[ <!-- --> ] -- .pull-right[ Всего возможно 15 сравнений. Если для каждого сравнения вероятность ошибки первого рода будет `\(\alpha_{per\ comparison} = 0.05\)`, то для группы из 15 сравнений --- ? ] -- Если предположить, что сравнения независимы (это не так), то `\(\alpha_{family\ wise} = 1 - (1 - 0.05) ^ {15} = 0.54\)`. Мы рискуем найти различия там где их нет с 54% вероятностью! Для зависимых сравнений вероятность будет немного меньше, но все равно значительно больше `\(0.05\)` --- ## Поправка Бонферрони --- очень жесткий способ коррекции. Если нужно много сравнений, можно снизить `\(\alpha _{per\ comparison}\)` до общепринятого уровня `$$\alpha _{per\ comparison} = \frac{\alpha _{family\ wise}}{n}$$` -- Например, если хотим зафиксировать `\(\alpha _{family\ wise} = 0.05\)` С поправкой Бонферрони `\(\alpha _{per\ comparison} = 0.05 / 15 = 0.003\)` Это очень жесткая поправка! Мы рискуем не найти достоверных различий, даже там, где они есть... Но есть выход. Вместо множества попарных сравнений можно использовать один тест --- дисперсионный анализ (analysis of variation, ANOVA). --- class: middle, center, inverse # Линейные модели с дискретными предикторами --- ## Для кодирования дискретных факторов в R используются две параметризации .pull-left[ __Параметризация индикаторных переменных__ Так же известна как - dummy coding - treatment parametrisation - reference cell model <br/> - В R обозначается __contr.treatment__. С ней вы уже знакомы. Используется по умолчанию в R. ] -- .pull-right[ __Параметризация эффектов__ <br/><br/> Так же известна как - effects coding - sum-to-zero parameterisation <br/><br/> - В R обозначается __contr.sum__. "Классическая" параметризация для дисперсионного анализа. Нужна, если хочется использовать т.наз. III тип сумм квадратов в многофакторном дисперсионном анализе со взаимодействием факторов. ] --- class: middle, center, inverse # Параметризация индикаторных переменных --- ## Переменные-индикаторы Фрагмент модельной матрицы: sp | spЛугКон <br/> `\(x_1\)` | spМалин <br/> `\(x_2\)` | spБелТряс <br/> `\(x_3\)` | spЛесКон <br/> `\(x_4\)` | spЛесЗав <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | | | | | ЛугКон | | | | | Малин | | | | | БелТряс | | | | | ЛесКон | | | | | ЛесЗав | | | | | --- ## Переменные-индикаторы Фрагмент модельной матрицы: sp | spЛугКон <br/> `\(x_1\)` | spМалин <br/> `\(x_2\)` | spБелТряс <br/> `\(x_3\)` | spЛесКон <br/> `\(x_4\)` | spЛесЗав <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | 0 | 0 | 0 | 0 | 0 ЛугКон | 1 | 0 | 0 | 0 | 0 Малин | 0 | 1 | 0 | 0 | 0 БелТряс | 0 | 0 | 1 | 0 | 0 ЛесКон | 0 | 0 | 0 | 1 | 0 ЛесЗав | 0 | 0 | 0 | 0 | 1 Переменных-индикаторов всегда на одну меньше, чем число уровней фактора. Уровень "Крапив" будет базовым: для его кодирования не нужна отдельная переменная. --- ## Коэффициенты модели в параметризации индикаторов .pull-left-60[ Фрагмент модельной матрицы: .small[ sp | spЛугКон <br/> `\(x_1\)` | spМалин <br/> `\(x_2\)` | spБелТряс <br/> `\(x_3\)` | spЛесКон <br/> `\(x_4\)` | spЛесЗав <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | 0 | 0 | 0 | 0 | 0 ЛугКон | 1 | 0 | 0 | 0 | 0 Малин | 0 | 1 | 0 | 0 | 0 БелТряс | 0 | 0 | 1 | 0 | 0 ЛесКон | 0 | 0 | 0 | 1 | 0 ЛесЗав | 0 | 0 | 0 | 0 | 1 ] ] .pull-right-40[ <!-- --> ] .pull-down[ `$$y _{i} = b_0 + b_1 x _{1i} + \ldots + b_5 x_{5i} + e_{i}$$` - `\(b_0\)` --- это среднее значение отклика для базового уровня фактора. - `\(b_1, ..., b_5\)` --- это отклонения от базового уровня для средних с другими уровнями фактора. ] --- ## Подбор модели в параметризации индикаторов ```r mod_treatment <- lm(len ~ sp, data = eggs) ``` ```r round(coef(mod_treatment), 2) ``` ``` (Intercept) spЛугКон spМалин spБелТряс spЛесКон spЛесЗав 21.12 1.17 1.44 1.77 1.96 1.99 ``` `$$\widehat{len}_i = 21.12 + 1.17 sp_{\text{ЛугКон}\ i} + 1.44 sp_{\text{Малин}\ i} + 1.77 sp_{\text{БелТряс}\ i} + \\ + 1.96 sp_{\text{ЛесКон}\ i} + 1.99 sp_{\text{ЛесЗав}\ i}$$` --- ## Уравнение модели в параметризации индикаторов `$$\widehat{len}_i = 21.12 + 1.17 sp_{\text{ЛугКон}\ i} + 1.44 sp_{\text{Малин}\ i} + 1.77 sp_{\text{БелТряс}\ i} + \\ + 1.96 sp_{\text{ЛесКон}\ i} + 1.99 sp_{\text{ЛесЗав}\ i}$$` Первый коэффициент --- средний размер яиц кукушек в гнездах крапивников (на базовом уровне): - `\(\widehat{len}_{\text{Крапив}} = 21.12\)` Другие коэффициенты --- разница размеров яиц кукушек в гнездах крапивников и других хозяев (отклонения от базового уровня). Если их прибавить к базовому уровню, то получим предсказанные значения для других видов: - `\(\widehat{len}_{\text{ЛугКон}\ i} = 21.12 + 1.17 sp_{\text{ЛугКон}\ i} = 22.29\)` - `\(\widehat{len}_{\text{Малин}\ i} = 21.12 + 1.44 sp_{\text{Малин}\ i} = 22.56\)` - `\(\widehat{len}_{\text{БелТряс}\ i} = 21.12 + 1.77 sp_{\text{БелТряс}\ i} = 22.89\)` - `\(\widehat{len}_{\text{ЛесКон}\ i} = 21.12 + 1.96 sp_{\text{ЛесКон}\ i}= 23.08\)` - `\(\widehat{len}_{\text{ЛесЗав}\ i} = 21.12 + 1.99 sp_{\text{ЛесЗав}\ i} = 23.11\)` --- class: middle, center, inverse # Параметризация эффектов --- ## Переменные-эффекты Фрагмент модельной матрицы: sp | sp1 <br/> `\(x_1\)` | sp2 <br/> `\(x_2\)` | sp3 <br/> `\(x_3\)` | sp4 <br/> `\(x_4\)` | sp5 <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | | | | | ЛугКон | | | | | Малин | | | | | БелТряс | | | | | ЛесКон | | | | | ЛесЗав | | | | | Переменных-эффектов всегда на одну меньше, чем число уровней фактора. Переменные закодированы при помощи -1, 0 и 1 (сумма кодов для возможных состояний одной переменной равна нулю). Для последней группы все переменные-эффекты будут равны `\(-1\)`. --- ## Переменные-эффекты Фрагмент модельной матрицы: sp | sp1 <br/> `\(x_1\)` | sp2 <br/> `\(x_2\)` | sp3 <br/> `\(x_3\)` | sp4 <br/> `\(x_4\)` | sp5 <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | 1 | 0 | 0 | 0 | 0 ЛугКон | 0 | 1 | 0 | 0 | 0 Малин | 0 | 0 | 1 | 0 | 0 БелТряс | 0 | 0 | 0 | 1 | 0 ЛесКон | 0 | 0 | 0 | 0 | 1 ЛесЗав | -1 | -1 | -1 | -1 | -1 Переменных-эффектов всегда на одну меньше, чем число уровней фактора. Переменные закодированы при помощи -1, 0 и 1 (сумма кодов для возможных состояний одной переменной равна нулю). Для последней группы все переменные-эффекты будут равны `\(-1\)`. --- ## Коэффициенты модели в параметризации эффектов .pull-left[ Фрагмент модельной матрицы: sp | sp1 <br/> `\(x_1\)` | sp2 <br/> `\(x_2\)` | sp3 <br/> `\(x_3\)` | sp4 <br/> `\(x_4\)` | sp5 <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | 1 | 0 | 0 | 0 | 0 ЛугКон | 0 | 1 | 0 | 0 | 0 Малин | 0 | 0 | 1 | 0 | 0 БелТряс | 0 | 0 | 0 | 1 | 0 ЛесКон | 0 | 0 | 0 | 0 | 1 ЛесЗав | -1 | -1 | -1 | -1 | -1 ] .pull-right[ <!-- --> ] .pull-down[ `$$y _{i} = b_0 + b_1 x _{1i} + \ldots + b_5 x_{5i} + e_{i}$$` - `\(b_0\)` --- это общее среднее значение отклика. - `\(b_1, ..., b_5\)` --- это отклонения от общего среднего для средних с другими уровнями фактора, кроме последнего. - для последнего уровня фактора отклонения от общего среднего --- это коэффициенты `\(b_1, ..., b_5\)`, взятые с противоположным знаком. ] --- ## Подбор модели в параметризации эффектов ```r mod_sum <- lm(len ~ sp, data = eggs, contrasts = list(sp = contr.sum)) ``` ```r round(coef(mod_sum), 2) ``` ``` (Intercept) sp1 sp2 sp3 sp4 sp5 22.51 -1.39 -0.22 0.05 0.38 0.57 ``` Коэффициенты моделей в разных параметризациях будут разными, но предсказания будут совершенно одинаковыми. `$$\widehat{len}_i = 22.51 -1.39 sp_{1\ i} -0.22 sp_{2\ i} + 0.05 sp_{3\ i} + 0.38 sp_{4\ i} + 0.57 sp_{5\ i}$$` --- ## Уравнение модели в параметризации эффектов `$$\widehat{len}_i = 22.51 -1.39 sp_{1\ i} -0.22 sp_{2\ i} + 0.05 sp_{3\ i} + 0.38 sp_{4\ i} + 0.57 sp_{5\ i}$$` Первый коэффициент --- средний размер яиц кукушек по всем данным: - `\(\overline{len} = 22.51\)` Другие коэффициенты — отличие размеров яиц в гнездах хозяев от общего среднего. Для всех уровней, кроме последнего — без изменения знака: - `\(\widehat{len}_{\text{Крапив}\ i} = 22.51 -1.39 sp_{1\ i} = 21.12\)` - `\(\widehat{len}_{\text{ЛугКон}\ i} = 22.51 -0.22 sp_{2\ i} = 22.29\)` - `\(\widehat{len}_{\text{Малин}\ i} = 22.51 + 0.05 sp_{3\ i} = 22.56\)` - `\(\widehat{len}_{\text{БелТряс}\ i} = 22.51 + 0.38 sp_{4\ i} = 22.89\)` - `\(\widehat{len}_{\text{ЛесКон}\ i} = 22.51 + 0.57 sp_{5\ i} = 23.08\)` Для последнего уровня фактора — с противоположным знаком, т.к. для него все переменные-эффекты `\(sp_1,..., sp_5 = -1\)`: - `\(\widehat{len}_{\text{ЛесЗав}\ i} = 22.51 -1.39 sp_{1\ i} -0.22 sp_{2\ i} + 0.05 sp_{3\ i} + 0.38 sp_{4\ i} + 0.57 sp_{5\ i} = 23.12\)` --- class: middle, center, inverse # t-тесты значимости коэффициентов --- ## t-тесты значимости коэффициентов не информативны в моделях с дискретными предикторами - Для модели __в параметризации индикаторов__ t-тесты коэффициентов показывают значимость отличий средних в группах от среднего на базовом уровне. - По значениям коэффициентов нельзя сказать влияет ли дискретный фактор целиком (исключение --- фактор с двумя градациями). ```r coef(summary(mod_treatment)) ``` ``` Estimate Std. Error t value Pr(>|t|) (Intercept) 21.120 0.2337 90.364 6.200e-108 spЛугКон 1.173 0.2699 4.348 3.007e-05 spМалин 1.436 0.3253 4.415 2.310e-05 spБелТряс 1.767 0.3305 5.345 4.699e-07 spЛесКон 1.960 0.3305 5.930 3.310e-08 spЛесЗав 1.994 0.3364 5.929 3.329e-08 ``` --- ## t-тесты значимости коэффициентов не информативны в моделях с дискретными предикторами - Для модели __в параметризации эффектов__ t-тесты коэффициентов показывают значимость отличий средних в группах от общего среднего -- это сравнение редко имеет смысл. ```r coef(summary(mod_sum)) ``` ``` Estimate Std. Error t value Pr(>|t|) (Intercept) 22.50842 0.09003 249.9974 5.090e-158 sp1 -1.38842 0.21101 -6.5800 1.492e-09 sp2 -0.21509 0.14229 -1.5117 1.334e-01 sp3 0.04783 0.20554 0.2327 8.164e-01 sp4 0.37824 0.21101 1.7926 7.569e-02 sp5 0.57158 0.21101 2.7088 7.794e-03 ``` --- class: segue-yellow # Дисперсионный анализ --- ## Общая изменчивость <!-- --> Общая изменчивость `\(SS_t\)` --- это сумма квадратов отклонений наблюдаемых значений `\(y_i\)` от общего среднего `\(\bar y\)` --- ## Факторная (межгрупповая) изменчивость <!-- --> Отклонения внутригрупповых средних от общего среднего в генеральной совокупности --- это эффект фактора `\(\alpha_j = \mu_j - \mu\)`, где `\(j = 1, 2, ..., p\)` --- это одна из `\(p\)` групп. Мы оцениваем эффект фактора по реальным данным `\(\bar{y}_j-\bar{y}\)` --- ## Структура общей изменчивости `$$SS_t = SS_x + SS_e$$` <!-- --> Общая изменчивость | Факторная изменчивость | Остаточная изменчивость ---- |----|---- ... | ... | ... `\(SS_{t}= \sum\sum{(\bar{y}-y_{ij})^2}\)` | `\(SS_{x}=\sum{n_j(\bar{y}_j-\bar{y})^2}\)` | `\(SS_{e}= \sum\sum{(\bar{y}_{j}-y_{ij})^2}\)` `\(df_{t} = n - 1\)` | `\(df_{x} = p - 1\)` | `\(df_{e} = n - p\)` --- ## От изменчивостей к дисперсиям `$$SS_t = SS_x + SS_e \qquad MS_t \ne MS_x + MS_e$$` <!-- --> Общая дисперсия | Факторная дисперсия | Остаточная дисперсия ---- |----|---- `\(MS_{t} = \frac {SS_{t}}{df_{t}}\)` | `\(MS_{x} = \frac {SS_{x}}{df_{x}}\)` | `\(MS_{e} = \frac{SS_{e}}{df_{e}}\)` `\(SS_{t}= \sum\sum{(\bar{y}-y_{ij})^2}\)` | `\(SS_{x}=\sum{n_j(\bar{y}_j-\bar{y})^2}\)` | `\(SS_{e}= \sum\sum{(\bar{y}_{j}-y_{ij})^2}\)` `\(df_{t} = n - 1\)` | `\(df_{x} = p - 1\)` | `\(df_{e} = n - p\)` ??? Они не зависят от числа наблюдений в выборке, в отличие от `\(SSx\)` и `\(SS_e\)` С их помощью можно проверить гипотезу о наличии связи между предиктором и откликом --- ## MSx и MSe <br/> помогают тестировать значимость фактора Если дисперсии остатков в группах равны и фактор имеет фиксированное число градаций: `\(E(MS_x) = \sigma^2 + \sum n_i \frac{(\mu_i - \mu)^2}{p - 1} = \sigma^2 + \sigma^2_x\)` `\(E(MS_e) = \sigma^2\)` -- Если зависимости нет, то `\(\mu_1 = \ldots = \mu_p\)` --- средние равны во всех `\(p\)` группах, и тогда `\(MS_x \sim MS_e\)`. -- - `\(H_0: \mu_1 = \ldots = \mu_p\)` --- средние во всех `\(p\)` группах равны. - `\(H_A: \exists\; i, j: \mu_i \ne \mu_j\)` --- __хотя бы одно__ среднее отличается от общего среднего. `$$F_{df_x, df_e} = \frac{MS _{x}}{MS_{e}}$$` --- ## Тестирование значимости фактора при помощи F-критерия - `\(H_0: \mu_1 = \ldots = \mu_p\)` --- средние во всех `\(p\)` группах равны. - `\(H_A: \exists\; i, j: \mu_i \ne \mu_j\)` --- __хотя бы одно__ среднее отличается от общего среднего. `$$F_{df_x, df_e} = \frac{MS _{x}}{MS_{e}}$$` В однофакторном дисперсионном анализе `\(df_{x} = p - 1\)` и `\(df_{e} = n - p\)`. <!-- --> --- ## Результаты дисперсионного анализа часто представляют в виде таблицы Источник изменчивости|SS|df|MS|F ----- | ----- | ----- | ----- | ----- Название фактора | `\(SS_{x}=\sum{n_j(\bar{y}_j-\bar{y})^2}\)` | `\(df _x = p - 1\)` | `\(MS _x = \frac{SS _x}{df _x}\)` | `\(F _{df _x df _e} = \frac{MS _x}{MS _e}\)` Случайная | `\(SS_{e}= \sum\sum{(\bar{y}_j-y_{ij})^2}\)` | `\(df _e = n - p\)` | `\(MS _e = \frac{SS _e}{df _e}\)` | Общая | `\(SS_{t}= \sum\sum{(\bar{y}-y_{ij})^2}\)` | `\(df _t = n - 1\)` | Минимальное описание результатов в тексте должно содержать `\(F _{df _x, df _e}\)` и `\(p\)`. --- ## Делаем дисперсионный анализ в R В R есть много функций для дисперсионного анализа. Мы рекомендуем `Anova()` (__с большой буквы__) из пакета `car`. Зачем? Эта функция умеет тестировать влияние факторов в определенном порядке. Когда факторов будет больше одного, это станет важно для результатов. ```r library(car) eggs_anova <- Anova(mod_treatment) eggs_anova ``` ``` Anova Table (Type II tests) Response: len Sum Sq Df F value Pr(>F) sp 42.8 5 10.4 0.000000029 *** Residuals 93.4 114 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Результаты дисперсионного анализа - Можно описать в тексте: Длина яиц кукушек в гнездах разных птиц-хозяев значимо различается <br/> ( `\(F _{ 5,114} = 10.45\)`, `\(p < 0.01\)`). - Можно представить в виде таблицы: Длина яиц кукушек значимо различалась в гнездах разных птиц-хозяев (Табл. 1). <table class="table" style="margin-left: auto; margin-right: auto;"> <caption>Табл. 1. Результаты дисперсионного анализа длины яиц кукушек в гнездах разных птиц-хозяев. SS --- суммы квадратов отклонений, df --- число степеней свободы, F --- значение F-критерия, P --- уровень значимости.</caption> <thead> <tr> <th style="text-align:left;"> Источник изменчивости </th> <th style="text-align:left;"> SS </th> <th style="text-align:right;"> df </th> <th style="text-align:left;"> F </th> <th style="text-align:left;"> P </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> Хозяин </td> <td style="text-align:left;"> 42.8 </td> <td style="text-align:right;"> 5 </td> <td style="text-align:left;"> 10.4 </td> <td style="text-align:left;"> &lt;0.01 </td> </tr> <tr> <td style="text-align:left;"> Остаточная </td> <td style="text-align:left;"> 93.4 </td> <td style="text-align:right;"> 114 </td> <td style="text-align:left;"> </td> <td style="text-align:left;"> </td> </tr> </tbody> </table> --- class: middle, center, inverse # Условия примененимости дисперсионного анализа --- ## Результатам тестов можно верить, если выполняются условия применимости Условия применимости дисперсионного анализа: - Случайность и независимость наблюдений внутри групп - Нормальное распределение остатков - Гомогенность дисперсий остатков - Отсутствие коллинеарности факторов (независимость групп) -- ### Другие ограничения: - Лучше работает, если размеры групп примерно одинаковы (т.наз. сбалансированный дисперсионный комплекс) - Устойчив к отклонениям от нормального распределения (при равных объемах групп или при больших выборках) --- ## Проверяем выполнение условий применимости ```r # Данные для графиков остатков mod_diag <- fortify(mod_treatment) ``` -- ### 1) График расстояния Кука ```r ggplot(mod_diag, aes(x = 1:nrow(mod_diag), y = .cooksd)) + geom_bar(stat = "identity") ``` <!-- --> -- - Выбросов нет --- ## 2) График остатков от предсказанных значений ```r ggplot(mod_diag, aes(x = .fitted, y = .stdresid)) + geom_jitter() ``` -- Но у нас один единственный дискретный предиктор, поэтому удобнее сразу боксплот остатков в зависимости от дискретного предиктора ### 3) Графики остатков от предикторов в модели и не в модели ```r ggplot(mod_diag, aes(x = sp, y = .stdresid)) + geom_boxplot() ``` <!-- --> -- - Дисперсии почти одинаковые. Может быть, в одной из групп чуть больше --- ## 4) Квантильный график остатков ```r library(car) qqPlot(mod_treatment, id = FALSE) ``` <!-- --> -- - Распределение остатков немного отличается от нормального --- class: middle, center, inverse # График предсказаний модели --- ## Данные для графика при помощи `predict()` ```r MyData <- data.frame(sp = factor(levels(eggs$sp), levels = levels(eggs$sp))) MyData <- data.frame( MyData, predict(mod_treatment, newdata = MyData, interval = "confidence")) MyData ``` ``` sp fit lwr upr 1 Крапив 21.12 20.66 21.58 2 ЛугКон 22.29 22.03 22.56 3 Малин 22.56 22.11 23.00 4 БелТряс 22.89 22.42 23.35 5 ЛесКон 23.08 22.62 23.54 6 ЛесЗав 23.11 22.64 23.59 ``` --- ## Задание 2 Создайте MyData вручную: - предсказанные значения - стандартные ошибки - верхнюю и нижнюю границы доверительных интервалов ```r MyData <- data.frame(sp = factor(levels(eggs$sp), levels = levels(eggs$sp))) X <- model.matrix() betas <- MyData$fit <- %*% MyData$se <- sqrt(diag(X %*% vcov(mod_treatment) %*% t(X))) t_crit <- qt(p = , df = nrow() - length(coef())) MyData$lwr <- MyData$ - * MyData$ MyData$upr <- MyData$ + * MyData$ ``` --- ## Решение 2 ```r MyData <- data.frame(sp = factor(levels(eggs$sp), levels = levels(eggs$sp))) X <- model.matrix(~sp, data = MyData) betas <- coef(mod_treatment) MyData$fit <- X %*% betas MyData$se <- sqrt(diag(X %*% vcov(mod_treatment) %*% t(X))) t_crit <- qt(p = 0.975, df = nrow(eggs) - length(coef(mod_treatment))) MyData$lwr <- MyData$fit - t_crit * MyData$se MyData$upr <- MyData$fit + t_crit * MyData$se MyData ``` ``` sp fit se lwr upr 1 Крапив 21.12 0.2337 20.66 21.58 2 ЛугКон 22.29 0.1349 22.03 22.56 3 Малин 22.56 0.2263 22.11 23.00 4 БелТряс 22.89 0.2337 22.42 23.35 5 ЛесКон 23.08 0.2337 22.62 23.54 6 ЛесЗав 23.11 0.2419 22.64 23.59 ``` --- ## Задание 3 Дополните код, чтобы получить столбчатый график с предсказаниями линейной модели. Заливкой покажите вид птиц-хозяев. Подпишите оси. Спрячьте легенду. ```r gg_bars <- ggplot(data = , aes(x = , y = )) + geom_col(aes(), width = 0.5) + geom_errorbar(aes(ymin = , ymax = ), width = 0.1) + labs( = "Вид хозяев", = "Длина яиц кукушек, мм") + scale_fill_brewer(name = "Вид \nхозяев", palette = "Dark2") + theme( = "none") gg_bars ``` <!-- --> --- ## Решение 3 Столбчатый график ```r gg_bars <- ggplot(data = MyData, aes(x = sp, y = fit)) + geom_col(aes(fill = sp), width = 0.5) + geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1) + labs(x = "Вид хозяев", y = "Длина яиц кукушек, мм") + scale_fill_brewer(name = "Вид \nхозяев", palette = "Dark2") + theme(legend.position = "none") gg_bars ``` <!-- --> -- Но какие именно из этих средних значений значимо различаются? --- class: segue-yellow # Пост хок тесты --- ## Как понять, какие именно группы различаются Дисперсионный анализ говорит нам только, есть ли влияние фактора, но не говорит, какие именно группы различаются. Коэффициенты линейной модели в `summary(mod_treatment)` содержат лишь часть ответа --- сравнение средних значених всех групп со средним на базовом уровне. Если нас интересуют другие возможные попарные сравнения, нужно сделать пост хок тест. --- ## Есть два способа понять, какие именно группы различаются -- .pull-left[ .center[__Линейные контрасты__ (linear contrasts)] - Гипотезы о межгрупповых различиях тестируются при помощи комбинаций из коэффициентов линейной модели. - Набор гипотез (и сравнений) должен быть определен заранее. - Делать можно вне зависимости от результатов дисперсионного анализа. Этот способ за рамками курса. ] -- .pull-right[ .center[__Post hoc__ тесты] - Сравниваются все возможные группы. - Нет четких заранее сформулированных гипотез. - Делать можно, только если влияние соответствующего фактора оказалось значимым. Этот способ мы обсудим. ] --- ## Разновидности пост хок тестов Тесты без поправки на число сравнений: - Наименьшая значимая разница Фишера (Fisher's Least Significant Difference) -- Тесты с поправкой для уровня значимости `\(\alpha\)`: - Поправка Бонферрони (Bonferroni correction) - Поправка Сидака (Sidak's correction) -- Тесты, основанные на распределении стьюдентизированного размаха: - Тест Тьюки (Tukey's Honest Significant Difference, HSD) - Тест Стьюдента-Ньюмена-Кьюлса (Student-Newman-Kewls test, SNK) - Тест Даннета (Dunnet's test) --- используется для сравнения с контрольной группой. -- Тесты, основанные на F-тестах: - Критерий Дункана (Dunkan's test) - Тест Шеффе (Scheffe's test) --- ## Наименьшая значимая разница Фишера <br/> Fisher's Least Significant Difference Используется t-критерий с `\(df = df_e = n - p\)`: <br/> `$$t = \frac{\bar{y}_i - \bar{y}_j}{\sqrt{MS_e \large(\frac{1}{n_i} + \frac{1}{n_j}\large)}}$$` -- - Подразумевается равенство дисперсий в сравниваемых группах - Не вносится поправка для уровня значимости, учитывающая множественность сравнений. (Считается, что тест "защищен" от ошибок I рода, т.к. выполняется после того, как в ANOVA была отвергнута гипотеза о равенстве всех внутригрупповых средних). -- __Осторожно!__ Этот тест слишком мягок, высока вероятность появления ошибок II рода (т.е. тест находит различия там, где их нет). --- ## После ANOVA часто приходится сравнивать несколько групп Фактор в дисперсионном анализе может задавать больше двух групп. (Например, фактор вид птицы-хозяина в нашем примере). -- На самом деле __t-распределение__ не годится для случая, когда приходится сравнивать больше, чем две группы одновременно. Вспомните, __t-распределение__ --- это распределение стандартизованной разницы <br/> средних значений __из двух выборок__, взятых из одной генеральной совокупности. -- <br/> Нужен способ описать __более сложное распределение --- для любого числа выборок__. --- ## Три выборки Представьте, что мы берем из одной и той же генеральной совокупности три выборки. Средние значения `\(\bar y_1\)`, `\(\bar y_2\)` и `\(\bar y_3\)` в каждой из этих выборок скорее всего окажутся разными и не будут похожи на генеральное среднее `\(\mu\)`. Как оценить, какой может быть эта разница? Нужно построить распределение. Но какое? -- <br/> 1.Возьмем `\(m\)` выборок из одной генеральной совокупности -- 2.Отсортируем выборочные средние: `\(\bar y_{1} \ge \bar y_2 \ge \ldots \ge \bar y_{m}\)` Это можно записать как `\(\bar y_{max} \ge \bar y_2 \ge \ldots \ge \bar y_{min}\)` -- 3.Вычислим разницу максимального и минимального средних `\(\bar y_{max} - \bar y_{min}\)` -- Если повторить 1--3 много раз, то получится распределение, которое показывает, чему может быть равна разница средних значений в выборках из одной генеральной совокупности. Такое распределение можно построить для любого числа выборок `\(m\)`. --- ## Распределение стьюдентизированного размаха <br/> Studentized range distribution Это распределение стандартизованной разницы минимального и максимального средних __для любого числа выборок__ из одной генеральной совокупности (форма зависит от `\(df\)` и от числа выборок `\(m\)`). .pull-left[ <!-- --> ] .pull-right[ Формула для случая равных дисперсий и разных объемов групп: `$$q = \frac{\bar{y}_{max} - \bar{y}_{min}}{\sqrt{s^2\frac{1}{2} \large(\frac{1}{n_i} + \frac{1}{n_j}\large)}}$$` ] --- ## Cтьюдентизированный t-критерий консервативнее обычного .pull-left[ .center[Обычный t-критерий] `$$t = \frac{\bar{y}_i - \bar{y}_j}{\sqrt{MS_e \large(\frac{1}{n_i} + \frac{1}{n_j}\large)}}$$` ] .pull-right[ .center[Стьюдентизированный t-критерий] `$$q = \frac{\bar{y}_i - \bar{y}_j}{\sqrt{MS_e \;\frac{1}{2} \; \large(\frac{1}{n_i} + \frac{1}{n_j}\large)}}$$` <!-- `$$q = \frac{\bar{y}_i - \bar{y}_j}{\sqrt{MS_e \;\tikzmarkin<2>{factor} \frac{1}{2} \tikzmarkend{factor} \; \large(\frac{1}{n_i} + \frac{1}{n_j}\large)}}$$`--> При этом `\(\bar{y}_i > \bar{y}_j\)`, т.е. вычитается из большего меньшее среднее. ] -- <br/> .center[Значение `\(q\)` будет в 1.414 раз больше, чем `\(t\)`. `$$q = \frac{t}{\sqrt{\; \frac{1}{2}}} = \sqrt{2} \cdot t = 1.414 \cdot t$$` ] --- ## Тест Тьюки (Tukey's Honest Significant Difference) Используется стьюдентизированный t-критерий с `\(df = df_e = n - p\)` и `\(m = p\)` (общее число групп): `$$q = \frac{\bar{y}_i - \bar{y}_j}{\sqrt{MS_e\frac{1}{2} \large(\frac{1}{n_i} + \frac{1}{n_j}\large)}}$$` Требуется равенство дисперсий. --- ## Пост хок тесты различаются по степени консервативности Если посмотреть на критические значения t при сравнении средних при `\(\alpha = 0.05\)` ($m = 4$ группы по 6 наблюдений, `\(df_e = 20\)`), становится понятно, что тест Тьюки --- разумный компромисс среди пост хок тестов. Тест | Критическое значение ---- | ----- Шеффе^a | 3.05 Бонферрони (4 группы) | 2.93 Тьюки (HSD)^b | 2.80 Бонферрони (3 группы) | 2.63 Даннет<sup>b<sup/> | 2.54 Дункан^{a, b} | 2.22 Фишер (LSD) | 2.09 ^a --- Значение `\(t\)` соответствующее `\(F\)`. ^b --- Для сопоставимости внесена поправка `\(\sqrt{2}\)`. --- ## Пост хок тест Тьюки в R - `glht()` --- "general linear hypotheses testing" - `linfct` --- аргумент, задающий гипотезу для тестирования - `mcp()` --- функция, чтобы задавать множественные сравнения (обычные пост хоки) - `sp` = "Tukey" --- тест Тьюки по фактору `sp` ```r library(multcomp) eggs_posthoc <- glht(mod_treatment, linfct = mcp(sp = "Tukey")) ``` --- ## Результаты попарных сравнений (тест Тьюки) Таблица результатов пост хок теста практически нечитабельна. ```r summary(eggs_posthoc) ``` ``` Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: lm(formula = len ~ sp, data = eggs) Linear Hypotheses: Estimate Std. Error t value Pr(>|t|) ЛугКон - Крапив == 0 1.1733 0.2699 4.35 <0.001 *** Малин - Крапив == 0 1.4362 0.3253 4.41 <0.001 *** БелТряс - Крапив == 0 1.7667 0.3305 5.34 <0.001 *** ЛесКон - Крапив == 0 1.9600 0.3305 5.93 <0.001 *** ЛесЗав - Крапив == 0 1.9943 0.3364 5.93 <0.001 *** Малин - ЛугКон == 0 0.2629 0.2635 1.00 0.915 БелТряс - ЛугКон == 0 0.5933 0.2699 2.20 0.241 ЛесКон - ЛугКон == 0 0.7867 0.2699 2.91 0.047 * ЛесЗав - ЛугКон == 0 0.8210 0.2770 2.96 0.041 * БелТряс - Малин == 0 0.3304 0.3253 1.02 0.909 ЛесКон - Малин == 0 0.5238 0.3253 1.61 0.587 ЛесЗав - Малин == 0 0.5580 0.3313 1.68 0.538 ЛесКон - БелТряс == 0 0.1933 0.3305 0.58 0.992 ЛесЗав - БелТряс == 0 0.2276 0.3364 0.68 0.984 ЛесЗав - ЛесКон == 0 0.0343 0.3364 0.10 1.000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) ``` --- ## Результаты пост хок теста Результаты пост хок теста можно привести в виде текста... - Размер яиц кукушек в гнездах крапивника значимо меньше, чем в гнездах лугового конька (тест Тьюки, `\(p < 0.01\)`). Размер яиц кукушек в гнездах лесной завирушки, белой трясогузки, малиновки и лесного конька не различается, но яйца кукушек в гнездах этих хозяев крупнее, чем в гнездах у лугового конька или крапивника (тест Тьюки, от `\(p < 0.01\)` до `\(0.05\)`). ...или построить график --- ## Можно привести результаты пост хок теста на столбчатом графике Значимо различающиеся группы обозначим разными буквами ```r gg_bars_coded <- gg_bars + geom_text(aes(y = 1.6, label = c("A", "B", "BC", "BC", "C", "C")), colour = "white", size = 7) gg_bars_coded ``` <!-- --> --- ## Take home messages - Дисперсионный анализ --- линейная модель с дискретными предикторами, существует в нескольких параметризациях, которые отличаются трактовками коэффициентов - При помощи дисперсионного анализа можно проверить гипотезу о равенстве средних значений в группах -- - Условия применимости дисперсионного анализа - Случайность и независимость групп и наблюдений внутри групп - Нормальное распределение в группах - Гомогенность дисперсий в группах -- - При множественных попарных сравнениях увеличивается вероятность ошибки первого рода, поэтому нужно вносить поправку для уровня значимости - Post hoc тесты --- это попарные сравнения после дисперсионного анализа, которые позволяют сказать, какие именно средние различаются --- ## Дополнительные ресурсы - Quinn, Keough, 2002, pp. 173--207 - Logan, 2010, pp. 254--282 - [Open Intro to Statistics](http://www.openintro.org/stat/), pp.236--246 - Sokal, Rohlf, 1995, pp. 179--260 - Zar, 2010, pp. 189-207