class: middle, left, inverse, title-slide .title[ # Дисперсионный анализ, часть 2 ] .subtitle[ ## Линейные модели… ] .author[ ### Марина Варфоломеева, Вадим Хайтов, Анастасия Лянгузова ] .date[ ### Осень 2024 ] --- ## Многофакторный дисперсионный анализ - Модель многофакторного дисперсионного анализа; - Взаимодействие факторов; - Несбалансированные данные, типы сумм квадратов; - Многофакторный дисперсионный анализ в R; - Дисперсионный анализ в матричном виде. ### Вы сможете - Проводить многофакторный дисперсионный анализ и интерпретировать его результаты с учётом взаимодействия факторов. --- class: middle, center, inverse # Данные --- ## Пример: Пингвины Зависит ли вес пингвина от его видовой принадлежности и пола? .pull-left[Измерения особей пингвинов из рода *Pygoscelis* лежат в датасете `penguins` в пакете `palmerpenguins`. Исходные данные были опубликованы в работе Gorman et al., 2014. Помимо веса и пола животных, датасет содержит информацию об острове, на котором пингвины проживали, и измерения клюва. В анализ мы возьмём только следующие переменные: Зависимая переменная: - `body_mass_g` --- вес в граммах. Факторы: - `species` --- вид пингвина; - `sex` --- пол пингвина.] .pull-right[  .small[ .pull-right[On [Pinterest](https://www.pinterest.com/pin/explore-the-natural-world--153544668520878653/)] ] ] --- ## Знакомимся с данными ``` r # Открываем данные library(palmerpenguins) peng <- as.data.frame(penguins[, c(1, 6, 7)]) str(peng) ``` ``` 'data.frame': 344 obs. of 3 variables: $ species : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ... $ body_mass_g: int 3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ... $ sex : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ... ``` ``` r # Переименовываем столбцы colnames(peng) <- c('sp', 'mass', 'sex') ``` --- ## Пропущенные значения ``` r colSums(is.na(peng)) ``` ``` sp mass sex 0 2 11 ``` Удаляем пропущенные значения (каких-то пингвинов не измерили :() ``` r pengs <- peng[complete.cases(peng), ] colSums(is.na(pengs)) ``` ``` sp mass sex 0 0 0 ``` --- ## Объемы выборок в группах ``` r table(pengs$sp, pengs$sex) ``` ``` female male Adelie 73 73 Chinstrap 34 34 Gentoo 58 61 ``` -- - Группы разного размера --- ## Посмотрим на график ``` r library(ggplot2) theme_set(theme_bw(base_size = 14)) gg_mass <- ggplot(data = pengs, aes(x = sp, y = mass, colour = sex)) + stat_summary(geom = 'pointrange', fun.data = mean_cl_normal) gg_mass ``` <!-- --> --- class: middle, center, inverse # Многофакторный дисперсионный анализ --- ## Многофакторный дисперсионный анализ Дисперсионный анализ становится многофакторным, если в модели используется несколько дискретных факторов. Дизайн такого рода анализа может быть разным. -- .pull-left[ ### Вложенные (иерархичные) модели В модели внутри какого-то фактора расположены дополнительные предикторы. Особенность таких моделей в том, что категории вложенных факторов в пределах уровней основного фактора различны. **Пример:** измерения обилия беспозвоночных животных на литорали в разные годы и месяцы. ] -- .pull-right[ ### Факторные модели Дизайн эксперимента предполагает, что разные факторы пересекаются друг с другом. Таким образом возможны любые комбинации между разными уровнями разных факторов. Наш случай! Факторов может быть несколько, мы будем разбираться на примере двухфакторного дисперсионного анализа. ] --- ## Какие бывают факторы Свойства | Фиксированные факторы | Случайные факторы ---- | ---- | ---- Уровни фактора | фиксированные, заранее определенные и потенциально воспроизводимые уровни | случайная выборка из всех возможных уровней Используются для тестирования гипотез | о средних значениях отклика между уровнями фактора <br/> `\(H _{0}: \mu _1 = \mu _2 = \ldots = \mu _i = \mu\)` | о дисперсии отклика между уровнями фактора <br/> `\(H _{0}: \sigma_{rand.fact.}^2 = 0\)` Выводы можно экстраполировать | только на уровни из анализа | на все возможные уровни Мы сейчас работает с фиксированными факторами. В случае, когда в модели присутстуют и фиксированные, и случайные факторы --- используют смешанные линейные модели, о которых мы поговорим позже. --- ## Взаимодействие факторов При добавлении новых факторов в модель появляется **взаимодействие** факторов. Взаимодействие факторов возникает, когда у одного фактора эффект разный в зависимости от уровней другого. Такое взаимодействие необходимо учитывать, в т.ч. при интепретации результатов. --- ## Что такое взаимодействие дискретных предикторов Взаимодействие факторов --- когда эффект фактора B разный в зависимости от уровней фактора A и наоборот. На каких рисунках есть взаимодействие факторов? (.small[Logan, 2010, fig.12.2]) .pull-left[  ] -- .pull-right[ - b, c - нет взаимодействия (эффект фактора B одинаковый для групп по фактору A, линии для разных групп по фактору B на графиках расположены параллельно) - a, d - есть взаимодействие (эффект фактора B разный для групп по фактору A, на графиках линии для разных групп по фактору B расположены под наклоном). ] --- ## Влияют ли главные эффекты и взаимодействие?  .small[Quinn, Keough, 2002, fig.9.3] -- - взаимодействие не значимо, и не мешает интерпретировать эффекты факторов. - фактор А влияет - фактор В влияет --- ## Влияют ли главные эффекты и взаимодействие?  .small[Quinn, Keough, 2002, fig.9.3] -- - взаимодействие значимо и мешает интерпретировать влияние факторов отдельно: - для В2 зависимая переменная возрастает с изменением уровня А - для В1 зависимая переменная возрастает только на А2, но не различается на А1 и А3 - __если смотреть на главные эффекты, можно сделать неправильные выводы (о факторе А)__: - фактор А влияет, группы А2 и А3 не отличаются - фактор В влияет, в группе В2 зависимая переменная больше, чем в В1 --- ## Влияют ли главные эффекты и взаимодействие?  .small[Quinn, Keough, 2002, fig.9.3] -- - взаимодействие значимо и мешает интерпретировать влияние факторов отдельно: - на уровне A2 меняется порядок различий уровней фактора B - __если смотреть на главные эффекты, можно сделать неправильные выводы__: - факторы А и В не влияют --- ## Взаимодействие факторов может маскировать главные эффекты .pull-left[  .small[Quinn, Keough, 2002, fig.9.3] ] .pull-right[ Если есть значимое взаимодействие, то - главные эффекты обсуждать не имеет смысла - пост хок тесты проводятся только для взаимодействия ] --- class: middle, center, inverse # Двухфакторный дисперсионный анализ <br/> в параметризации индикаторов --- ## Переменные-индикаторы В нашем примере отклик --- вес пингвина, и два дискретных фактора: - `sex` --- 2 уровня (базовый `female`), для кодирования нужна одна переменная. ``` r contr.treatment(levels(pengs$sex)) ``` ``` male female 0 male 1 ``` - `sp` --- 3 уровня (базовый `Adelie`), для кодирования нужно две переменных. ``` r contr.treatment(levels(pengs$sp)) ``` ``` Chinstrap Gentoo Adelie 0 0 Chinstrap 1 0 Gentoo 0 1 ``` --- ## Переменные-индикаторы Дополнительные переменные понадобятся, чтобы учесть взаимодействие факторов. Фрагмент модельной матрицы: sex <br/> | sp <br/> | sexsp <br/> `\(x_1\)` | spChinstrap <br/> `\(x_2\)` | spGentoo <br/> `\(x_3\)`| sexmale:Chinstrap <br/> `\(x_4\)` | sexmale:Gentoo <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- | ---- female | Adelie | 0 | 0 | 0 | 0 | 0 male | Adelie | 1 | 0 | 0 | 0 | 0 female | Chinstrap | 0 | 1 | 0 | 0 | 0 male | Chinstrap | 1 | 1 | 0 | 1 | 0 female | Gentoo | 0 | 0 | 1 | 0 | 0 male | Gentoo | 1 | 0 | 1 | 0 | 1 --- ## Уравнение линейной модели в параметризации индикаторов `$$y _{i} = b _0 + b _1 x _{1i} + b _2 x _{2i} + b _3 x _{3i} + b _4 x _{4i} + b _5 x _{5i}+ e _{i}$$` - `\(b_0\)` --- значение отклика для самок вида Adelie (на базовом уровне обоих факторов). Отклонения относительно базового уровня обоих факторов: - `\(b_1\)` --- для самцов Adelie; - `\(b_2\)` и `\(b_3\)` --- для самок Chinstrap и Gentoo, соответственно; - `\(b_4\)` и `\(b_5\)` --- для самцов Chinstrap и Gentoo, соответственно. --- ## Подбираем линейную модель <br/> в параметризации индикаторов (contr.treatment) ``` r mod_treat_pengs <- lm(mass ~ sex * sp, data = pengs) mod_treat_pengs ``` ``` Call: lm(formula = mass ~ sex * sp, data = pengs) Coefficients: (Intercept) sexmale spChinstrap spGentoo 3369 675 158 1311 sexmale:spChinstrap sexmale:spGentoo -263 130 ``` Общее уравнение модели `$$\begin{aligned}\widehat{mass} _{i} = 3369 + 675 sex_{male,i} + 158 sp_{Chinstrap\,i} + 1311 sp_{Gentoo\,i} + \\ - 263 sex_{male}\ sp_{Chinstrap\,i} + 130 sex_{male}\ sp_{Gentoo\,i} \end{aligned}$$` --- class: middle, center, inverse # Двухфакторный дисперсионный анализ в параметризации эффектов --- ## Переменные-эффекты В нашем примере отклик --- вес пингвина, и два дискретных фактора: - `sex` --- 2 уровня (базовый `female`), для кодирования нужна одна переменная. ``` r contr.sum(levels(pengs$sex)) ``` ``` [,1] female 1 male -1 ``` - `sp` --- 3 уровня (базовый `Adelie`), для кодирования нужно две переменных. ``` r contr.sum(levels(pengs$sp)) ``` ``` [,1] [,2] Adelie 1 0 Chinstrap 0 1 Gentoo -1 -1 ``` --- ## Переменные-эффекты Дополнительные переменные понадобятся, чтобы учесть взаимодействие факторов. Фрагмент модельной матрицы: sex <br/> | sp <br/> | sex1 <br/> `\(x_1\)` | sp1 <br/> `\(x_2\)` | sp2 <br/> `\(x_3\)` | sex1:sp1 <br/> `\(x_4\)` | sex1:sp2 <br/> `\(x_5\)` :---- | :----: | :----: | :----: | :----: | :----: | :----: female | Adelie | 1 | 1 | 0 | 1 | 0 male | Adelie | -1 | 1 | 0 | -1 | 0 female | Chinstrap | 1 | 0 | 1 | 0 | 1 male | Chinstrap | -1 | 0 | 1 | 0 | -1 female | Gentoo | 1 | -1 | -1 | 1 | 1 male | Gentoo | -1 | -1 | -1 | -1 | -1 --- ## Уравнение линейной модели в параметризации эффектов `$$y _{i} = b _0 + b _1 x _{1i} + b _2 x _{2i} + b _3 x _{3i} + b _4 x _{4i} + b _5 x _{5i}+ e _{i}$$` - `\(b_0\)` --- среднее значение отклика по всем данным. Отклонения от общего среднего значений отклика: - `\(b_1\)` --- в зависимости от пола (фактор `sex`); - `\(b_2\)` и `\(b_3\)` --- в зависимости от вида (фактор `sp`); - `\(b_4\)` и `\(b_5\)` --- для пола в зависимости от вида (взаимодействие). --- ## Подбираем линейную модель <br/> в параметризации эффектов (contr.sum) ``` r mod_sum_pengs <- lm(mass ~ sex * sp, data = pengs, contrasts = list(sp = 'contr.sum', sex = 'contr.sum')) coef(mod_sum_pengs) ``` ``` (Intercept) sex1 sp1 sp2 sex1:sp1 sex1:sp2 4173.85 -315.25 -467.68 -440.76 -22.08 109.37 ``` Общее уравнение модели `$$\begin{aligned}\widehat{mass}_i = 417.85 - 315.25 sex_{female\,i} - 467.68 sp_{Adelie\,i} - 440.76 sp_{Chinstrap\,i} - \\ - 22.08 sex_{female\,i}sp_{Adelie\,i} + 109.37 sex_{female\,i}sp_{Chinstrap\,i} \end{aligned}$$` --- class: middle, center, inverse # Диагностика линейной модели --- ## Диагностика линейной модели Нужно проверить, выполняются ли условия применимости <br/> для модели в нужной параметризации Данные для анализа остатков ``` r mod_treat_peng_diag <- fortify(mod_treat_pengs) # функция из пакета ggplot2 head(mod_treat_peng_diag, 2) ``` ``` mass sex sp .hat .sigma .cooksd .fitted .resid 1 3750 male Adelie 0.0137 309.4 0.002112 4043 -293.5 2 3800 female Adelie 0.0137 308.9 0.004558 3369 431.2 .stdresid 1 -0.9552 2 1.4032 ``` --- ## График расстояния Кука ``` r ggplot(mod_treat_peng_diag, aes(x = 1:nrow(mod_treat_peng_diag), y = .cooksd)) + geom_bar(stat = 'identity') ``` <!-- --> -- - Влиятельных наблюдений нет. --- ## График остатков от предсказанных значений ``` r gg_resid <- ggplot(data = mod_treat_peng_diag, aes(x = .fitted, y = .stdresid)) + geom_point() + geom_hline(yintercept = 0) gg_resid ``` <!-- --> -- - Влиятельных наблюдений нет (все в пределах 3 SD). - Несколько наблюдений не совсем обычны (2SD < .stdresid < 3SD), но поскольку расстояние Кука для них небольшое, это нестрашно. --- ## График зависимости остатков от предикторов в модели ``` r ggplot(data = mod_treat_peng_diag, aes(x = sex, y = .stdresid, colour = sp)) + geom_boxplot() + geom_hline(yintercept = 0) ``` <!-- --> Удобнее смотреть на боксплот. Нет гетерогенности дисперсии, всё хорошо! --- ## Квантильный график остатков ``` r library(car) qqPlot(mod_treat_pengs, id = FALSE) # функция из пакета car ``` <!-- --> -- - Отклонений от нормального распределения нет. --- class: middle, center, inverse # Несбалансированные данные, типы сумм квадратов --- ## Несбалансированные данные - когда численности в группах по факторам различаются .pull-left[ Например так, | | A1 | A2 | A3 | |----|----|----|----| | B1 | 5 | 5 | 5 | | B2 | 5 | 4 | 5 | ] .pull-right[ или так, | | A1 | A2 | A3 | |----|----|----|----| | B1 | 3 | 8 | 4 | | B2 | 4 | 7 | 4 | ] --- ## Проблемы из-за несбалансированности данных - Оценки средних в разных группах с разным уровнем точности (Underwood 1997) - ANOVA менее устойчив к отклонениям от условий применимости (особенно от __гомогенности__ дисперсий) при разных размерах групп (Quinn Keough 2002, section 8.3) - Проблемы с расчетом мощности. Если `\(\sigma _{\epsilon}^2 > 0\)` и размеры выборок разные, то `\(MS _{x} \over MS _{e}\)` не следует F-распределению (Searle et al. 1992). <br/> -- Старайтесь _планировать_ группы равной численности! Но если не получилось --- не страшно: - Для фиксированных эффектов неравные размеры --- проблема при нарушении условий применимости только, если значения доверительной вероятности _p_ близки к выбранному критическому уровню значимости `\(\alpha\)`. --- ## Суммы квадратов в многофакторном дисперсионном анализе со взаимодействием ### Если данные сбалансированы, то ... - взаимодействие и эффекты факторов независимы (в любой параметризации), - все суммы квадратов и соответствующие тесты можно посчитать в одном анализе, - результат не зависит от порядка включения факторов в модель. ### Если данные несбалансированы, то ... - суммы квадратов для факторов не равны общей сумме квадратов, - для вычислений используется регрессионный подход (несколько сравнений вложенных моделей), - результат анализа может зависеть от порядка включения факторов в модель. --- ## Порядок тестирования значимости предикторов в дисперсионном анализе "Типы сумм квадратов" | I тип | II тип | III тип ---- | ---- | ---- | ---- Название | Последовательный | Без учета взаимодействий высоких порядков | Иерархический --- ## Порядок тестирования значимости предикторов <br/> в дисперсионном анализе .small[ "Типы сумм квадратов" | I тип | II тип | III тип ---- | ---- | ---- | ---- Название | Последовательный | Без учета взаимодействий высоких порядков | Иерархический Порядок расчета SS | SS(A) <br/> SS(B&#124;A) <br/> SS(AB&#124;B, A) | SS(A&#124;B) <br/> SS(B&#124;A) <br/> SS(AB&#124;B, A) | SS(A&#124;B, AB) <br/> SS(B&#124;A, AB) <br/> SS(AB&#124;B, A) Величина эффекта зависит от выборки в группе | Да | Да | Нет Результат зависит от порядка включения факторов в модель | Да | Нет | Нет Параметризация | Любая | Любая | Только параметризация эффектов Команда R | aov(), anova() | Anova() (пакет car) | Anova() (пакет car) __Осторожно!__ Тестируя предикторы в разном порядке, вы тестируете разные гипотезы! ] --- ## Если несбалансированные данные, выберите подходящий порядок тестирования гипотез <!-- - SSe и SSab всегда рассчитываются одинаково, вне зависимости от порядка тестирования гипотез и от сбалансированности данных --> <!-- - SSa, SSb --- есть три способа расчета (суммы квадратов I, II и III типа, терминология пришла из SAS) в зависимости от порядка тестирования значимости факторов --> ### Если данные сбалансированы, то ... - При использовании любого типа сумм квадратов результаты расчетов будут одинаковы. ### Если данные несбалансированы, то ... - Результаты зависят от выбранного типа сумм квадратов (т.к. он определяет, какие гипотезы при этом тестируются). <br/> Для несбалансированных данных иногда рекомендуют __суммы квадратов III типа__ если есть взаимодействие факторов (Maxwell & Delaney 1990, Milliken, Johnson 1984, Searle 1993, Yandell 1997, Glantz, Slinker 2000). Но при этом __нарушается принцип маргинальности__, поэтому некоторые статистики не любят тех, кто так делает... --- class: middle, center, inverse # Многофакторный дисперсионный анализ в R --- ## Дисперсионный анализ со II типом сумм квадратов При таком способе, сначала тестируется взаимодействие, затем отдельные факторы в модели без взаимодействия. ``` r mod_treat_pengs <- lm(mass ~ sp * sex, data = pengs) library(car) Anova(mod_treat_pengs, type = "II") ``` ``` Anova Table (Type II tests) Response: mass Sum Sq Df F value Pr(>F) sp 143401584 2 749.02 <2e-16 *** sex 37090262 1 387.46 <2e-16 *** sp:sex 1676557 2 8.76 0.0002 *** Residuals 31302628 327 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Дисперсионный анализ c III типом сумм квадратов Опишем процедуру на тот случай, если вдруг вам понадобится воспроизвести в R дисперсионный анализ с III типом сумм квадратов. При этом способе вначале тестируют взаимодействие, когда все другие факторы есть в модели. Затем тестируют факторы, когда все другие факторы и взаимодействие есть в модели. -- __Внимание: при использовании III типа сумм квадратов, нужно обязательно указывать тип контрастов для факторов__ (`contrasts=list(фактор_1 = contr.sum, фактор_2=contr.sum)`). -- ``` r mod_sum_pengs <- lm(mass ~ sp * sex, data = pengs, contrasts = list(sp = 'contr.sum', sex = 'contr.sum')) Anova(mod_sum_pengs, type = "III") ``` ``` Anova Table (Type III tests) Response: mass Sum Sq Df F value Pr(>F) (Intercept) 5232595969 1 54661.83 <2e-16 *** sp 143001222 2 746.92 <2e-16 *** sex 29851220 1 311.84 <2e-16 *** sp:sex 1676557 2 8.76 0.0002 *** Residuals 31302628 327 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Почему для расчета III типа сумм квадратов обязательно использовать параметризацию эффектов ? Для расчета III типа сумм квадратов нужно иметь возможность удалить из модели влияние предиктора, и одновременно оставить в ней взаимодействие (т.е. предикторы и взаимодействие были независимы друг от друга). -- __В параметризации индикаторных переменных предикторы и взаимодействие коллинеарны__, т.е. суммы квадратов III типа будут рассчитаны неправильно. ``` r vif(mod_treat_pengs) ``` ``` GVIF Df GVIF^(1/(2*Df)) sp 4.066 2 1.420 sex 2.281 1 1.510 sp:sex 6.688 2 1.608 ``` -- __В параметризации эффектов переменных предикторы и взаимодействие независимы__, значит получатся верные суммы квадратов III типа. ``` r vif(mod_sum_pengs) ``` ``` GVIF Df GVIF^(1/(2*Df)) sp 1.000 2 1.000 sex 1.109 1 1.053 sp:sex 1.109 2 1.026 ``` --- class: middle, center, inverse # Пост хок тест для взаимодействия факторов --- ## Пост хок тесты в многофакторном дисперсионном анализе - Поскольку взаимодействие достоверно, факторы отдельно можно не тестировать. Проведем пост хок тест по взаимодействию, чтобы выяснить, какие именно группы различаются - Если бы взаимодействие было недостоверно, мы бы провели пост хок тест по тем факторам, влияние которых было бы достоверно. Как? См. предыдущую презентацию. --- ## Пост хок тест для взаимодействия факторов Пост хок тест для взаимодействия факторов делается легче всего "обходным путём": 1. Создаем переменную-взаимодействие; 2. Подбираем модель без свободного члена; 3. Делаем пост хок тест для этой модели. --- ## Задание 1 Дополните этот код, чтобы посчитать пост хок тест Тьюки по взаимодействию факторов ``` r # Создаем переменную-взаимодействие pengs$sp_sex <- interaction(pengs$sex, pengs$sp) # Подбираем линейную модель без свободного члена fit_inter <- lm() # Делаем пост хок тест для этой модели library() dat_tukey <- glht(, linfct = mcp( = )) summary(dat_tukey) ``` --- ## Решение ``` r # Создаем переменную-взаимодействие pengs$sex_sp <- interaction(pengs$sex, pengs$sp) # Подбираем линейную модель без свободного члена fit_inter <- lm(mass ~ sex_sp - 1, data = pengs) # Делаем пост хок тест для этой модели library(multcomp) dat_tukey <- glht(fit_inter, linfct = mcp(sex_sp = 'Tukey')) summary(dat_tukey) ``` --- ## Результаты пост хок теста в виде таблицы почти нечитабельны ``` Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: lm(formula = mass ~ sex_sp - 1, data = pengs) Linear Hypotheses: Estimate Std. Error t value Pr(>|t|) male.Adelie - female.Adelie == 0 674.7 51.2 13.17 <0.001 *** female.Chinstrap - female.Adelie == 0 158.4 64.2 2.47 0.14 male.Chinstrap - female.Adelie == 0 570.1 64.2 8.88 <0.001 *** female.Gentoo - female.Adelie == 0 1310.9 54.4 24.09 <0.001 *** male.Gentoo - female.Adelie == 0 2116.0 53.7 39.43 <0.001 *** female.Chinstrap - male.Adelie == 0 -516.3 64.2 -8.04 <0.001 *** male.Chinstrap - male.Adelie == 0 -104.5 64.2 -1.63 0.58 female.Gentoo - male.Adelie == 0 636.2 54.4 11.69 <0.001 *** male.Gentoo - male.Adelie == 0 1441.3 53.7 26.85 <0.001 *** male.Chinstrap - female.Chinstrap == 0 411.8 75.0 5.49 <0.001 *** female.Gentoo - female.Chinstrap == 0 1152.5 66.8 17.25 <0.001 *** male.Gentoo - female.Chinstrap == 0 1957.6 66.2 29.56 <0.001 *** female.Gentoo - male.Chinstrap == 0 740.8 66.8 11.08 <0.001 *** male.Gentoo - male.Chinstrap == 0 1545.9 66.2 23.35 <0.001 *** male.Gentoo - female.Gentoo == 0 805.1 56.7 14.19 <0.001 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) ``` --- ## Данные для графика при помощи `predict()` У нас два дискретных фактора, поэтому вначале используем `expand.grid()` ``` r MyData_pengs <- expand.grid(sex = levels(pengs$sex), sp = levels(pengs$sp)) MyData_pengs <- data.frame( MyData_pengs, predict(mod_treat_pengs, newdata = MyData_pengs, interval = 'confidence') ) MyData_pengs ``` ``` sex sp fit lwr upr 1 female Adelie 3369 3298 3440 2 male Adelie 4043 3972 4115 3 female Chinstrap 3527 3423 3632 4 male Chinstrap 3939 3835 4043 5 female Gentoo 4680 4600 4760 6 male Gentoo 5485 5407 5563 ``` --- ## Задание 2 Создайте MyData вручную для модели в обычной параметризации: - предсказанные значения - стандартные ошибки - верхнюю и нижнюю границы доверительных интервалов ``` r MyData_pengs <- expand.grid(sex = levels(pengs$sex), sp = levels(pengs$sp)) X_pengs <- model.matrix(~ , data = ) betas_pengs <- coef() MyData_pengs$fit <- MyData_pengs$se <- (( %*% vcov() %*% t())) MyData_pengs$lwr <- MyData_pengs$ - 2 * MyData_pengs$ MyData_pengs$upr <- MyData_pengs$ + 2 * MyData_pengs$ MyData_pengs ``` ``` sex sp fit se lwr upr 1 female Adelie 3369 36.21 3296 3441 2 male Adelie 4043 36.21 3971 4116 3 female Chinstrap 3527 53.06 3421 3633 4 male Chinstrap 3939 53.06 3833 4045 5 female Gentoo 4680 40.63 4598 4761 6 male Gentoo 5485 39.61 5406 5564 ``` --- ## Решение: ``` r MyData_pengs <- expand.grid(sex = levels(pengs$sex), sp = levels(pengs$sp)) X_pengs <- model.matrix(~ sp * sex, data = MyData_pengs) betas_pengs <- coef(mod_treat_pengs) MyData_pengs$fit <- X_pengs %*% betas_pengs MyData_pengs$se <- sqrt(diag(X_pengs %*% vcov(mod_treat_pengs) %*% t(X_pengs))) MyData_pengs$lwr <- MyData_pengs$fit - 2 * MyData_pengs$se MyData_pengs$upr <- MyData_pengs$fit + 2 * MyData_pengs$se MyData_pengs ``` ``` sex sp fit se lwr upr 1 female Adelie 3369 36.21 3296 3441 2 male Adelie 4043 36.21 3971 4116 3 female Chinstrap 3527 53.06 3421 3633 4 male Chinstrap 3939 53.06 3833 4045 5 female Gentoo 4680 40.63 4598 4761 6 male Gentoo 5485 39.61 5406 5564 ``` --- ## Задание 3 Постройте график результатов, на котором будут изображены предсказанные средние значения видового богатства в зависимости от тритмента и времени экспозиции. ``` r pos <- position_dodge(width = 0.2) gg_linep <- ggplot(data = , aes()) + geom_ (position = pos) + geom_ (aes(group = ), position = pos) + geom_ (position = pos, width = 0.1) gg_linep ``` <!-- --> --- ## График результатов: Линии с точками ``` r pos <- position_dodge(width = 0.2) gg_linep <- ggplot(data = MyData_pengs, aes(x = sp, y = fit, ymin = lwr, ymax = upr, colour = sex)) + geom_point(position = pos) + geom_line(aes(group = sex), position = pos) + geom_errorbar(position = pos, width = 0.1) gg_linep ``` <!-- --> --- ## Приводим график в приличный вид ``` r gg_final <- gg_linep + labs(x = 'Вид', y = 'Масса тела, г') + scale_colour_manual(name = '', labels = c('Самки', 'Самцы'), values = c('#FF0091', '#0077FF')) gg_final ``` <!-- --> --- ## Take home messages - Многофакторный дисперсионный анализ позволяет оценить взаимодействие факторов. Если оно значимо, то лучше воздержаться от интерпретации их индивидуальных эффектов -- - Если численности групп равны, получаются одинаковые результаты вне зависимости от порядка тестирования значимости факторов -- - В случае, если численности групп неравны (несбалансированные данные), есть несколько способов тестирования значимости факторов (I, II, III типы сумм квадратов) --- ## Дополнительные ресурсы - Quinn, Keough, 2002, pp. 221-250 - Logan, 2010, pp. 313-359 - Sokal, Rohlf, 1995, pp. 321-362 - Zar, 2010, pp. 246-266