class: middle, left, inverse, title-slide .title[ # Линейные модели с дискретными предикторами ] .subtitle[ ## Линейные модели… ] .author[ ### Марина Варфоломеева, Вадим Хайтов, Анастасия Лянгузова ] .date[ ### Осень 2024 ] --- # Линейные модели с дискретными предикторами (дисперсионный анализ) ## Вы сможете - Объяснить, в чем опасность множественных сравнений, и как с ними можно бороться; - Рассказать, как в дисперсионном анализе моделируются значения зависимой переменной; - Интерпретировать и описать результаты, записанные в таблице дисперсионного анализа; - Перечислить и проверить условия применимости дисперсионного анализа; - Провести множественные попарные сравнения при помощи 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 # Параметризация индикаторных переменных --- ## Переменные-индикаторы Основная идея --- присвоить уникальные значения различным значениям фактора, используя __дихотомию__ 0 и 1. Записать это можно в форме модельной матрицы. sp | spКрапив <br/> `\(x_1\)` | spЛугКон <br/> `\(x_2\)` | spМалин <br/> `\(x_3\)` | spБелТряс <br/> `\(x_4\)` | spЛесКон <br/> `\(x_5\)` | spЛесЗав <br/> `\(x_6\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | 1 | 0 | 0 | 0 | 0 | 0 ЛугКон | 0 | 1 | 0 | 0 | 0 | 0 Малин | 0 | 0 | 1 | 0 | 0 | 0 БелТряс | 0 | 0 | 0 | 1 | 0 | 0 ЛесКон | 0 | 0 | 0 | 0 | 1 | 0 ЛесЗав | 0 | 0 | 0 | 0 | 0 | 1 Такая форма записи избыточна! Можем принять первый фактор за референсное значение. --- ## Переменные-индикаторы Модельную матрицу параметризации записывают следующим образом: 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\)` ---- | ---- | ---- | ---- | ---- | ---- Крапив | | | | | ЛугКон | | | | | Малин | | | | | БелТряс | | | | | ЛесКон | | | | | ЛесЗав | | | | | --- ## Переменные-эффекты Числа в пределах модельной матрицы в случае параметризации эффектов будут образовывать __трихотомию__ из 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 Единица кодирует принадлежность к конкретной категории фактора, 0 --- принадлежность к другому уровню, начиная от первого (Крапивник в нашем случае) и заканчивая предпоследним уровнем. У последней группы все переменные-эффекты будут равны `\(-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[ <!-- --> ] <br /> `$$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: middle, center, inverse # Дисперсионный анализ --- ## Общая изменчивость <!-- --> Общая изменчивость `\(SS_t\)` --- это сумма квадратов отклонений наблюдаемых значений `\(y_i\)` от общего среднего `\(\bar y\)` --- ## Отклонения от общего среднего .pull-left[ **Межгрупповые отклонения** --- отклонения внутригрупповых средних от общего среднего ("эффекты" групп) --- факторная изменчивость. <!-- --> ] .pull-right[ **Внутригрупповые отклонения** --- отклонения наблюдаемых значений от внутригрупповых средних --- случайная изменчивость. <br> <!-- --> ] Отклонения внутригрупповых средних от общего среднего в генеральной совокупности --- это эффект фактора `\(\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$$` <!-- --> Общая дисперсия | Факторная дисперсия | Остаточная дисперсия ---- |----|---- `\(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\)` `\(MS_{t} = \frac {SS_{t}}{df_{t}}\)` | `\(MS_{x} = \frac {SS_{x}}{df_{x}}\)` | `\(MS_{e} = \frac{SS_{e}}{df_{e}}\)` ??? Они не зависят от числа наблюдений в выборке, в отличие от `\(SSx\)` и `\(SS_e\)` С их помощью можно проверить гипотезу о наличии связи между предиктором и откликом --- ## `\(MS_x\)` и `\(MS_e\)` <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;"> <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: middle, center, inverse # Пост хок тесты --- ## Как понять, какие именно группы различаются Дисперсионный анализ говорит нам только, есть ли влияние фактора, но не говорит, какие именно группы различаются. Коэффициенты линейной модели в `summary(mod_treatment)` содержат лишь часть ответа --- сравнение средних значених всех групп со средним на базовом уровне. Если нас интересуют другие возможные попарные сравнения, нужно сделать пост хок тест. --- ## Есть два способа понять, какие именно группы различаются .pull-left[ .center[__Линейные контрасты__ (linear contrasts)] - Гипотезы о межгрупповых различиях тестируются при помощи комбинаций из коэффициентов линейной модели. - Набор гипотез (и сравнений) должен быть определен заранее. - Делать можно вне зависимости от результатов дисперсионного анализа. ] -- .pull-right[ .center[__Post hoc__ тесты] - Сравниваются все возможные группы. - Нет четких заранее сформулированных гипотез. - Делать можно, только если влияние соответствующего фактора оказалось значимым. ] --- ## Разновидности пост хок тестов Тесты без поправки на число сравнений: - Наименьшая значимая разница Фишера (Fisher's Least Significant Difference) -- Тесты с поправкой для уровня значимости `\(\alpha\)`: - Поправка Бонферрони (Bonferroni correction) - Поправка Сидака (Šidák'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 часто приходится сравнивать несколько групп Фактор в дисперсионном анализе может задавать больше двух групп. (Например, фактор вид птицы-хозяина в нашем примере). -- На самом деле __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{MS_{e}\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)}}$$` Требуется равенство дисперсий. --- ## Тест Даннета (Dunnett’s test) Используется, когда хотим сравнить контрольную группу с экспериментальными. Лучше подходит для тех случаев, когда мы не хотим сравнить все группы со всеми, а нам нужны только определённые. Основан также на стьюдентизированном t-критерии. q-статистика считается так же, как в случае теста Тьюки; по-другому расчитывается _критическое_ значение q, с которым сравниваем наблюдаемое. --- ## Пост хок тесты различаются по степени консервативности Если посмотреть на критические значения t при сравнении средних при `\(\alpha = 0.05\)` `\((m = 4\)` группы по 6 наблюдений, `\(df_e = 20\)`), становится понятно, что тест Тьюки --- разумный компромисс среди пост хок тестов. .pull-left[ <!-- --> ] .pull-right[ Тест | Критическое значение ---- | ----- Шеффе<sup>a</sup> | 3.05 Бонферрони (4 группы) | 2.93 Тьюки (HSD)<sup>b</sup> | 2.80 Бонферрони (3 группы) | 2.63 Даннет<sup>b<sup/> | 2.54 Дункан<sup>a, b</sup> | 2.22 Фишер (LSD) | 2.09 ] <sup>a</sup> --- Значение `\(t\)` соответствующее `\(F\)`. <sup>b</sup> --- Для сопоставимости внесена поправка `\(\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 ``` <!-- --> --- class: middle, center, inverse ## Линейные контрасты --- ## Применение линейных контрастов Когда нам нужны только какие-то определённые группы сравнения, а не все группы со всеми сразу. Модель сравнения записывается в форме модельной матрицы. Для начала разберёмся со сравнением двуx групп в модели, заданной параметризацией индикаторов. Допустим, мы хотим сравнить лугового конька с белой трясогузкой. Возьмём лугового конька в качестве референсной группы. В таком случае можно оценить контраст для неё как: `\(L_1 = 1 * b_{ЛугКон} = 1 * b_{ЛугКон} + 0 * b_{БелТряс}\)`. Данное значение будет идентично интерсепту в R. Контраст для второй группы сравнения можно записать как разница между средним значением в ней и референсным значением: `\(L_2 = 1 * b_{БелТряс} - 1 * b_{ЛугКон}\)`. Таким образом, получаем матрицу следующего вида: <div class="math"> \[ \begin{array}{cc} & \begin{vmatrix} 1 & 0 \\ 1 & -1 \end{vmatrix}\end{array} \] </div> --- ## Линейные контрасты Так же, как в случае записи модельной матрицы, такая полная запись несколько излишня, поскольку нас интересует только та часть матрицы, которая соответствует __сравнению__ между двумя группами. Т.е. БелТряс --- ЛугКон = <div class="math"> \[ \begin{array}{cc} & \begin{vmatrix} 1 & -1 \end{vmatrix}\end{array} \] </div> __НО!__ В нашем датасете гораздо больше групп сравнения, и нас может интересовать большее количество групп. Например, хотим ещё сравнить размер лугового конька с лесной завирушкой. Уравнение для данного случая можно записать следующим образом: `$$\textit{БелТряс} - \textit{ЛугКон} = 0$$` `$$\textit{ЛесЗав} - \textit{ЛугКон} = 0$$` --- ## Запись линейных контрастов текстом в R ``` r eggs_contrtext <- glht(mod_treatment, linfct = mcp(sp = c("БелТряс - ЛугКон = 0", "ЛесЗав - ЛугКон = 0"))) summary(eggs_contrtext) ``` ``` Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: lm(formula = len ~ sp, data = eggs) Linear Hypotheses: Estimate Std. Error t value Pr(>|t|) БелТряс - ЛугКон == 0 0.593 0.270 2.20 0.0580 . ЛесЗав - ЛугКон == 0 0.821 0.277 2.96 0.0073 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) ``` --- ## Запись матрицы контрастов с несколькими группами сравнения В нашиx данныx несколько видов птиц, сначала нужно записать уравнение с учётом каждого значения фактора. Уравнение для длины яиц в гнездах белой трясогузки: `$$Y_{БелТряс} = 1 * b_0 + 0 * b_{ЛугКон} + 0 * b_{Малин} + \\ + 1 * b_{БелТряс} + 0 * b_{ЛесКон} + 0 * b_{ЛесЗав}$$` в гнездах лугового конька: `$$Y_{ЛугКон} = 1 * b_0 + 1 * b_{ЛугКон} + 0 * b_{Малин} + \\ + 0 * b_{БелТряс} + 0 * b_{ЛесКон} + 0 * b_{ЛесЗав}$$` в гнездах лесной завирушки: `$$Y_{ЛесЗав} = 1 * b_0 + 0 * b_{ЛугКон} + 0 * b_{Малин} + \\ + 0 * b_{БелТряс} + 0 * b_{ЛесКон} + 1 * b_{ЛесЗав}$$` Для сравнения групп так же вычитаем члены одного уравнения из другого. --- ## Запись матрицы контрастов с несколькими группами сравнения Например, для сравнения белой трясогузки и лугового конька. `$$БелТряс - ЛугКон = 0 * b_0 - 1 * b_{ЛугКон} + 0 * b_{Малин} + \\ + 1 * b_{БелТряс} + 0 * b_{ЛесКон} + 0 * b_{ЛесЗав}$$` Лесной завирушки и лугового конька: `$$ЛесЗав - ЛугКон = 0 * b_0 - 1 * b_{ЛугКон} + 0 * b_{Малин} + \\ + 0 * b_{БелТряс} + 0 * b_{ЛесКон} + 1 * b_{ЛесЗав}$$` --- ## Запись матрицы контрастов в R В случае когда у нас, например, две группы сравнения: `\(\textit{БелТряс} - \textit{ЛугКон}\)` и `\(\textit{ЛесЗав} - \textit{ЛугКон}\)`, матрица контрастов будет выглядеть на наших данных следующим образом. <div class="math"> \[ \begin{array}{cc} & \begin{vmatrix} 0 & -1 & 0 & 1 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & 1 \end{vmatrix}\end{array} \] </div> Код будет выглядеть следующим образом: ``` r contr <- rbind("БелТряс - ЛугКон" = c(0, -1, 0, 1, 0, 0), "ЛесЗав - ЛугКон" = c(0, -1, 0, 0, 0, 1)) eggs_contrmat <- glht(mod_treatment, linfct = contr) summary(eggs_contrmat) ``` ``` Simultaneous Tests for General Linear Hypotheses Fit: lm(formula = len ~ sp, data = eggs) Linear Hypotheses: Estimate Std. Error t value Pr(>|t|) БелТряс - ЛугКон == 0 0.593 0.270 2.20 0.0580 . ЛесЗав - ЛугКон == 0 0.821 0.277 2.96 0.0073 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) ``` --- ## 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