class: middle, left, inverse, title-slide .title[ # Дисперсионный анализ, часть 2 ] .subtitle[ ## Линейные модели… ] .author[ ### Марина Варфоломеева, Вадим Хайтов ] .date[ ### Осень 2022 ] --- ## Многофакторный дисперсионный анализ - Модель многофакторного дисперсионного анализа - Взаимодействие факторов - Несбалансированные данные, типы сумм квадратов - Многофакторный дисперсионный анализ в R - Дисперсионный анализ в матричном виде ### Вы сможете - Проводить многофакторный дисперсионный анализ и интерпретировать его результаты с учетом взаимодействия факторов --- class: middle, center, inverse # Данные --- ## Пример: Удобрение и беспозвоночные Влияет ли добавление азотных и фосфорных удобрений на беспозвоночных? Небольшие искуственные субстраты экспонировали в течение разного времени в верхней части сублиторали (Hall et al., 2000). Зависимая переменная: - `richness` --- Число видов Факторы: - `time` --- срок экспозиции (2, 4 и 6 месяцев) - `treat` --- удобрения (добавляли или нет) Планировали сделать 5 повторностей для каждого сочетания факторов .small[Данные:Quinn, Keough, 2002] --- ## Знакомимся с данными ```r fert <- read.csv(file='data/hall.csv') str(fert) ``` ``` 'data.frame': 29 obs. of 3 variables: $ TREAT : chr "control" "control" "control" "control" ... $ TIME : int 2 2 2 2 2 4 4 4 4 4 ... $ RICHNESS: int 5 7 5 7 5 20 18 20 18 17 ... ``` ```r # Для удобства названия переменных маленькими буквами colnames(fert) <- tolower(colnames(fert)) # Факторы делаем факторами fert$treat <- factor(fert$treat) fert$time <- factor(fert$time) ``` --- ## Пропущенные значения ```r colSums(is.na(fert)) ``` ``` treat time richness 0 0 0 ``` -- - Нет пропущенных значений --- ## Объемы выборок в группах ```r table(fert$time, fert$treat) ``` ``` control nutrient 2 5 5 4 5 5 6 4 5 ``` -- - Группы разного размера --- ## Посмотрим на график ```r library(ggplot2) theme_set(theme_bw(base_size = 14)) gg_rich <- ggplot(data = fert, aes(x = time, y = richness, colour = treat)) + stat_summary(geom = 'pointrange', fun.data = mean_cl_normal) gg_rich ``` <!-- --> -- - Вполне возможно, здесь есть гетерогенность дисперсий. Позже мы это проверим. --- ## Преобразовываем данные Зависимая переменная `richness` -- это счетная величина. Она подчиняется распределению Пуассона (и чем больше ее среднее значение, тем больше дисперсия). Правильно было бы воспользоваться обобщенными линейными моделями с Пуассоновским распределением ошибок вместо нормального. Но сейчас мы с вами попробуем действовать грубее (пока еще не разобрались, как это делать правильно). Давайте мы попробуем преобразовать зависимую переменную, чтобы ее распределение стало больше походить на нормальное. Это может помочь, а может и нет. ```r fert$log_rich <- log10(fert$richness + 1) ``` --- class: middle, center, inverse # Многофакторный дисперсионный анализ --- ## Многофакторный дисперсионный анализ Дисперсионный анализ становится многофакторным, если в модели используется несколько дискретных факторов. В таком анализе появляется взаимодействие факторов. Взаимодействие факторов возникает, когда у одного фактора эффект разный в зависимости от уровней другого. Разберемся с этим на схемах. --- ## Что такое взаимодействие дискретных предикторов Взаимодействие факторов - когда эффект фактора 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/> в параметризации индикаторов --- ## Переменные-индикаторы В нашем примере отклик --- видовое богатство, и два дискретных фактора: -- - `treat` --- 2 уровня (базовый `control`), для кодирования нужна одна переменная. ```r contr.treatment(levels(fert$treat)) ``` ``` nutrient control 0 nutrient 1 ``` -- - `time` --- 3 уровня (базовый `2`), для кодирования нужно две переменных. ```r contr.treatment(levels(fert$time)) ``` ``` 4 6 2 0 0 4 1 0 6 0 1 ``` --- ## Переменные-индикаторы Дополнительные переменные понадобятся, чтобы учесть взаимодействие факторов. Фрагмент модельной матрицы: treat <br/> | time <br/> | treatnutrient <br/> `\(x_1\)` | time4 <br/> `\(x_2\)` | time6 <br/> `\(x_3\)`| treatnutrient:time4 <br/> `\(x_4\)` | treatnutrient:time6 <br/> `\(x_5\)` ---- | ---- | ---- | ---- | ---- | ---- | ---- control | 2 | 0 | 0 | 0 | 0 | 0 nutrient | 2 | 1 | 0 | 0 | 0 | 0 control | 4 | 0 | 1 | 0 | 0 | 0 nutrient | 4 | 1 | 1 | 0 | 1 | 0 control | 6 | 0 | 0 | 1 | 0 | 0 nutrient | 6 | 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\)` --- значение отклика для контроля через 2 месяца (на базовом уровне обоих факторов) Отклонения относительно базового уровня обоих факторов: - `\(b_1\)` --- для удобренных площадок - `\(b_2\)` и `\(b_3\)` --- для площадок с экспозицией 4 или 6 мес - `\(b_4\)` и `\(b_5\)` --- для удобренных площадок с экспозицией 4 или 6 мес --- ## Подбираем линейную модель <br/> в параметризации индикаторов (contr.treatment) ```r mod_treatment <- lm(log_rich ~ treat * time, data = fert) mod_treatment ``` ``` Call: lm(formula = log_rich ~ treat * time, data = fert) Coefficients: (Intercept) treatnutrient time4 time6 0.8281 0.0504 0.4633 0.6031 treatnutrient:time4 treatnutrient:time6 0.1381 0.0462 ``` Общее уравнение модели `$$\begin{aligned}\widehat{log\_rich} _{i} = 0.828 + 0.05 treat_{nutrient\,i} + 0.463 time_{4\,i} + 0.603 time_{6\,i} + \\ + 0.138 treat_{nutrient}\ time_{4\,i} + 0.046 treat_{nutrient}\ time_{6\,i} \end{aligned}$$` --- class: middle, center, inverse # Двухфакторный дисперсионный анализ в параметризации эффектов --- ## Переменные-эффекты В нашем примере отклик --- видовое богатство, и два дискретных фактора: -- - `treat` --- 2 уровня (базовый `control`), для кодирования нужна одна переменная. ```r contr.sum(levels(fert$treat)) ``` ``` [,1] control 1 nutrient -1 ``` -- - `time` --- 3 уровня (базовый `2`), для кодирования нужно две переменных. ```r contr.sum(levels(fert$time)) ``` ``` [,1] [,2] 2 1 0 4 0 1 6 -1 -1 ``` --- ## Переменные-эффекты Дополнительные переменные понадобятся, чтобы учесть взаимодействие факторов. Фрагмент модельной матрицы: treat | time | treat1 <br/> `\(x_1\)` | time1 <br/> `\(x_2\)` | time2 <br/> `\(x_3\)` | treat1:time1 <br/> `\(x_4\)` | treat1:time2 <br/> `\(x_5\)` :---- | :----: | :----: | :----: | :----: | :----: | :----: control | 2 | 1 | 1 | 0 | 1 | 0 nutrient | 2 | -1 | 1 | 0 | -1 | 0 control | 4 | 1 | 0 | 1 | 0 | 1 nutrient | 4 | -1 | 0 | 1 | 0 | -1 control | 6 | 1 | -1 | -1 | -1 | -1 nutrient | 6 | -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\)` --- в зависимости от тритмента (фактор `treat`) - `\(b_2\)` и `\(b_3\)` --- в зависимости от экспозиции (фактор `time`) - `\(b_4\)` и `\(b_5\)` --- для тритментов в зависимости от экспозиции (взаимодействие) --- ## Подбираем линейную модель <br/> в параметризации эффектов (contr.sum) ```r mod_sum <- lm(log_rich ~ treat * time, data = fert, contrasts = list(treat = 'contr.sum', time = 'contr.sum')) mod_sum ``` ``` Call: lm(formula = log_rich ~ treat * time, data = fert, contrasts = list(treat = "contr.sum", time = "contr.sum")) Coefficients: (Intercept) treat1 time1 time2 treat1:time1 treat1:time2 1.2395 -0.0559 -0.3862 0.1462 0.0307 -0.0383 ``` Общее уравнение модели `$$\begin{aligned}\widehat{log\_rich}_i = 1.24 -0.056 treat_{1\,i} -0.386 time_{1\,i} + 0.146 time_{2\,i} + \\ + 0.031 treat_{1\,i}time_{1\,i} -0.038 treat_{1\,i}time_{2\,i} \end{aligned}$$` --- class: middle, center, inverse # Диагностика линейной модели --- ## Диагностика линейной модели Нужно проверить, выполняются ли условия применимости <br/> для модели в нужной параметризации Данные для анализа остатков ```r mod_treatment_diag <- fortify(mod_treatment) # функция из пакета ggplot2 head(mod_treatment_diag, 2) ``` ``` log_rich treat time .hat .sigma .cooksd .fitted 1 0.7782 control 2 0.2 0.05671 0.04050 0.8281 2 0.9031 control 2 0.2 0.05512 0.09112 0.8281 .resid .stdresid 1 -0.04998 -0.9859 2 0.07496 1.4788 ``` --- ## График расстояния Кука ```r ggplot(mod_treatment_diag, aes(x = 1:nrow(mod_treatment_diag), y = .cooksd)) + geom_bar(stat = 'identity') ``` <!-- --> -- - Влиятельных наблюдений нет. --- ## График остатков от предсказанных значений ```r gg_resid <- ggplot(data = mod_treatment_diag, aes(x = .fitted, y = .stdresid)) + geom_point() + geom_hline(yintercept = 0) gg_resid ``` <!-- --> -- - Влиятельных наблюдений нет (все в пределах 3 SD). --- ## График зависимости остатков от предикторов в модели .pull-left[ ```r ggplot(data = mod_treatment_diag, aes(x = treat, y = .stdresid, colour = time)) + geom_boxplot() + geom_hline(yintercept = 0) ``` <!-- --> Удобнее смотреть на боксплот. ] -- .pull-right[ - Видна гетерогенность дисперсии. В данном случае это не страшно, т.к. дисперсионный анализ устойчив к ситуации, когда в одной из групп разброс меньше, чем в других (особенно, если данные не слишком несбалансированные) (Underwood, 1997, McGuinness, 2002) ] --- ## Квантильный график остатков ```r library(car) qqPlot(mod_treatment, id = FALSE) # функция из пакета car ``` <!-- --> -- - Отклонений от нормального распределения нет. --- class: segue-yellow # Несбалансированные данные, типы сумм квадратов --- ## Несбалансированные данные - когда численности в группах по факторам различаются .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_treatment <- lm(log_rich ~ treat * time, data = fert) library(car) Anova(mod_treatment, type = 'II') ``` ``` Anova Table (Type II tests) Response: log_rich Sum Sq Df F value Pr(>F) treat 0.091 1 28.42 0.000021 *** time 2.216 2 344.88 < 2e-16 *** treat:time 0.025 2 3.84 0.036 * Residuals 0.074 23 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Дисперсионный анализ c III типом сумм квадратов .small[ Опишем процедуру на тот случай, если вдруг вам понадобится воспроизвести в R дисперсионный анализ с III типом сумм квадратов. При этом способе вначале тестируют взаимодействие, когда все другие факторы есть в модели. Затем тестируют факторы, когда все другие факторы и взаимодействие есть в модели. ] -- __Внимание: при использовании III типа сумм квадратов, нужно обязательно указывать тип контрастов для факторов__ (`contrasts=list(фактор_1 = contr.sum, фактор_2=contr.sum)`). -- ```r mod_sum <- lm(log_rich ~ treat * time, data = fert, contrasts = list(treat = contr.sum, time = contr.sum)) Anova(mod_sum, type = 3) ``` ``` Anova Table (Type III tests) Response: log_rich Sum Sq Df F value Pr(>F) (Intercept) 44.2 1 13776.11 < 2e-16 *** treat 0.1 1 28.04 0.000022 *** time 2.2 2 344.75 < 2e-16 *** treat:time 0.0 2 3.84 0.036 * Residuals 0.1 23 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Почему для расчета III типа сумм квадратов обязательно использовать параметризацию эффектов ? Для расчета III типа сумм квадратов нужно иметь возможность удалить из модели влияние предиктора, и одновременно оставить в ней взаимодействие (т.е. предикторы и взаимодействие были независимы друг от друга). -- __В параметризации индикаторных переменных предикторы и взаимодействие коллинеарны__, т.е. суммы квадратов III типа будут рассчитаны неправильно. ```r vif(mod_treatment) ``` ``` GVIF Df GVIF^(1/(2*Df)) treat 2.897 1 1.702 time 4.345 2 1.444 treat:time 8.517 2 1.708 ``` -- __В параметризации эффектов переменных предикторы и взаимодействие независимы__, значит получатся верные суммы квадратов III типа. ```r vif(mod_sum) ``` ``` GVIF Df GVIF^(1/(2*Df)) treat 1.006 1 1.003 time 1.009 2 1.002 treat:time 1.009 2 1.002 ``` --- class: segue-yellow # Пост хок тест для взаимодействия факторов --- ## Пост хок тесты в многофакторном дисперсионном анализе - Поскольку взаимодействие достоверно, факторы отдельно можно не тестировать. Проведем пост хок тест по взаимодействию, чтобы выяснить, какие именно группы различаются - Если бы взаимодействие было недостоверно, мы бы провели пост хок тест по тем факторам, влияние которых было бы достоверно. Как? См. предыдущую презентацию. --- ## Пост хок тест для взаимодействия факторов Пост хок тест для взаимодействия факторов делается легче всего "обходным путем" 1. Создаем переменную-взаимодействие 2. Подбираем модель без свободного члена 3. Делаем пост хок тест для этой модели --- ## Задание 1 Дополните этот код, чтобы посчитать пост хок тест Тьюки по взаимодействию факторов ```r # Создаем переменную-взаимодействие fert$treat_time <- interaction(fert$treat, fert$time) # Подбираем линейную модель от этой переменной без свободного члена fit_inter <- lm() # Делаем пост хок тест для этой модели library(multcomp) dat_tukey <- glht(, linfct = mcp( = 'Tukey')) summary() ``` --- ## Решение ```r # Создаем переменную-взаимодействие fert$treat_time <- interaction(fert$treat, fert$time) # Подбираем линейную модель без свободного члена fit_inter <- lm(log_rich ~ treat_time - 1, data = fert) # Делаем пост хок тест для этой модели library(multcomp) dat_tukey <- glht(fit_inter, linfct = mcp(treat_time = 'Tukey')) summary(dat_tukey) ``` --- ## Результаты пост хок теста в виде таблицы почти нечитабельны ``` Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: lm(formula = log_rich ~ treat_time - 1, data = fert) Linear Hypotheses: Estimate Std. Error t value Pr(>|t|) nutrient.2 - control.2 == 0 0.0504 0.0358 1.41 0.723 control.4 - control.2 == 0 0.4633 0.0358 12.93 <0.001 *** nutrient.4 - control.2 == 0 0.6519 0.0358 18.19 <0.001 *** control.6 - control.2 == 0 0.6031 0.0380 15.86 <0.001 *** nutrient.6 - control.2 == 0 0.6997 0.0358 19.52 <0.001 *** control.4 - nutrient.2 == 0 0.4129 0.0358 11.52 <0.001 *** nutrient.4 - nutrient.2 == 0 0.6015 0.0358 16.78 <0.001 *** control.6 - nutrient.2 == 0 0.5527 0.0380 14.54 <0.001 *** nutrient.6 - nutrient.2 == 0 0.6493 0.0358 18.11 <0.001 *** nutrient.4 - control.4 == 0 0.1885 0.0358 5.26 <0.001 *** control.6 - control.4 == 0 0.1398 0.0380 3.68 0.014 * nutrient.6 - control.4 == 0 0.2364 0.0358 6.59 <0.001 *** control.6 - nutrient.4 == 0 -0.0488 0.0380 -1.28 0.791 nutrient.6 - nutrient.4 == 0 0.0478 0.0358 1.33 0.763 nutrient.6 - control.6 == 0 0.0966 0.0380 2.54 0.153 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) ``` --- ## Данные для графика при помощи `predict()` У нас два дискретных фактора, поэтому вначале используем `expand.grid()` ```r MyData <- expand.grid(treat = levels(fert$treat), time = levels(fert$time)) MyData <- data.frame( MyData, predict(mod_treatment, newdata = MyData, interval = 'confidence') ) # Обратная трансформация (не забываем про единичку, которую прибавляли) MyData$richness <- 10^MyData$fit - 1 MyData$LWR <- 10^MyData$lwr - 1 MyData$UPR <- 10^MyData$upr - 1 MyData ``` ``` treat time fit lwr upr richness LWR UPR 1 control 2 0.8281 0.7757 0.8806 5.732 4.966 6.596 2 nutrient 2 0.8785 0.8261 0.9310 6.560 5.700 7.530 3 control 4 1.2914 1.2390 1.3439 18.563 16.339 21.074 4 nutrient 4 1.4800 1.4276 1.5324 29.199 25.764 33.073 5 control 6 1.4312 1.3726 1.4898 25.991 22.583 29.891 6 nutrient 6 1.5278 1.4754 1.5802 32.713 28.879 37.040 ``` --- ## Задание 2 Создайте MyData вручную для модели в обычной параметризации: - предсказанные значения - стандартные ошибки - верхнюю и нижнюю границы доверительных интервалов ```r MyData <- expand.grid(treat = levels(fert$treat), time = levels()) X <- model.matrix(~ , data = ) betas <- coef() MyData$fit <- MyData$se <- (X %*% vcov(mod_treatment) %*% t(X)) MyData$lwr <- MyData$ - 2 * MyData$upr <- MyData$ + 2 * # Обратная трансформация MyData$richness <- MyData$LWR <- MyData$UPR <- MyData ``` ``` treat time fit se lwr upr richness LWR UPR 1 control 2 0.8281 0.02535 0.7774 0.8788 5.732 4.990 6.565 2 nutrient 2 0.8785 0.02535 0.8278 0.9292 6.560 5.727 7.496 3 control 4 1.2914 0.02535 1.2408 1.3421 18.563 16.408 20.985 4 nutrient 4 1.4800 0.02535 1.4293 1.5307 29.199 25.872 32.937 5 control 6 1.4312 0.02834 1.3745 1.4879 25.991 22.689 29.753 6 nutrient 6 1.5278 0.02535 1.4771 1.5785 32.713 28.999 36.887 ``` --- ## Решение: ```r MyData <- expand.grid(treat = levels(fert$treat), time = levels(fert$time)) X <- model.matrix(~ treat * time, data = MyData) betas <- coef(mod_treatment) MyData$fit <- X %*% betas MyData$se <- sqrt(diag(X %*% vcov(mod_treatment) %*% t(X))) MyData$lwr <- MyData$fit - 2 * MyData$se MyData$upr <- MyData$fit + 2 * MyData$se # Обратная трансформация MyData$richness <- 10^MyData$fit - 1 MyData$LWR <- 10^MyData$lwr - 1 MyData$UPR <- 10^MyData$upr - 1 MyData ``` ``` treat time fit se lwr upr richness LWR UPR 1 control 2 0.8281 0.02535 0.7774 0.8788 5.732 4.990 6.565 2 nutrient 2 0.8785 0.02535 0.8278 0.9292 6.560 5.727 7.496 3 control 4 1.2914 0.02535 1.2408 1.3421 18.563 16.408 20.985 4 nutrient 4 1.4800 0.02535 1.4293 1.5307 29.199 25.872 32.937 5 control 6 1.4312 0.02834 1.3745 1.4879 25.991 22.689 29.753 6 nutrient 6 1.5278 0.02535 1.4771 1.5785 32.713 28.999 36.887 ``` --- ## Задание 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, aes(x = time, y = richness, ymin = LWR, ymax = UPR, colour = treat)) + geom_point(position = pos) + geom_line(aes(group = treat), position = pos) + geom_errorbar(position = pos, width = 0.1) gg_linep ``` <!-- --> --- ## Приводим график в приличный вид ```r gg_final <- gg_linep + labs(x = 'Экспозиция', y = 'Число видов') + scale_colour_brewer(name = '', palette = 'Dark2', labels = c('Контроль', 'Эксперимент')) 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