class: middle, left, inverse, title-slide .title[ # Дисперсионный анализ, часть 2 ] .subtitle[ ## Линейные модели… ] .author[ ### Марина Варфоломеева, Вадим Хайтов, Анастасия Лянгузова ] .date[ ### Осень 2025 ] --- ## Многофакторный дисперсионный анализ - Модель многофакторного дисперсионного анализа; - Взаимодействие факторов; - Несбалансированные данные, типы сумм квадратов; - Многофакторный дисперсионный анализ в R; - Дисперсионный анализ в матричном виде. ### Вы сможете - Проводить многофакторный дисперсионный анализ и интерпретировать его результаты с учётом взаимодействия факторов. --- class: middle, center, inverse # One Way ANOVA Update: Линейные контрасты --- ## Открываем данные ``` 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("ЛесЗав", "ЛугКон", "БелТряс", "Малин", "ЛесКон", "Крапив") ``` --- ## Как понять, какие именно группы различаются Дисперсионный анализ говорит нам только, есть ли влияние фактора, но не говорит, какие именно группы различаются. Коэффициенты линейной модели в `summary(mod_treatment)` содержат лишь часть ответа --- сравнение средних значених всех групп со средним на базовом уровне. Если нас интересуют другие возможные попарные сравнения, нужно сделать пост хок тест. --- ## Есть два способа понять, какие именно группы различаются -- .pull-left[ .center[__Post hoc__ тесты] - Сравниваются все возможные группы. - Нет четких заранее сформулированных гипотез. - Делать можно, только если влияние соответствующего фактора оказалось значимым. Этот способ мы уже обсуждали на предыдущей лекции. ] -- .pull-right[ .center[__Линейные контрасты__ (linear contrasts)] - Гипотезы о межгрупповых различиях тестируются при помощи комбинаций из коэффициентов линейной модели. - Набор гипотез (и сравнений) **должен быть определен заранее**. - Делать можно вне зависимости от результатов дисперсионного анализа. Этот способ за рамками курса. Однако в некоторых отраслях это важно, мы обсудим лишь общую суть линейных контрастов. ] --- ## Применение линейных контрастов Когда нам нужны только какие-то определённые группы сравнения, а не все группы со всеми сразу. Модель сравнения записывается в форме модельной матрицы. Для начала разберёмся со сравнением дву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 mod_treatment <- lm(len ~ sp, data = eggs) library(multcomp) 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.263 0.263 1.00 0.53 ЛесЗав - ЛугКон == 0 -1.173 0.270 -4.35 0.00006 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) ``` --- class: middle, center, inverse # One Way ANOVA Update: параметризация модели --- ## Для кодирования дискретных факторов в 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 - В R обозначается __contr.sum__. "Классическая" параметризация для дисперсионного анализа. Нужна, если изначально предполагается, что факторы взаимодействуют и взаимодействия мы удалять не собираемся. В этой ситуации придется использовать III тип сумм квадратов. - SS(A | B, AB) для фактора A. - SS(B | A, AB) для фактора B - SS(AB | B, A) для взаимодействия ] --- 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Крапив 23.11 -0.82 -0.23 -0.56 -0.03 -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{Крапив}} = 23.11\)` Другие коэффициенты --- разница размеров яиц кукушек в гнездах крапивников и других хозяев (отклонения от базового уровня). Если их прибавить к базовому уровню, то получим предсказанные значения для других видов: - `\(\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.88\)` - `\(\widehat{len}_{\text{БелТряс}\ i} = 21.12 + 1.77 sp_{\text{БелТряс}\ i} = 22.55\)` - `\(\widehat{len}_{\text{ЛесКон}\ i} = 21.12 + 1.96 sp_{\text{ЛесКон}\ i}= 23.08\)` - `\(\widehat{len}_{\text{ЛесЗав}\ i} = 21.12 + 1.99 sp_{\text{ЛесЗав}\ i} = 21.12\)` --- 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 0.61 -0.22 0.38 0.05 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} = 23.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.89\)` - `\(\widehat{len}_{\text{БелТряс}\ i} = 22.51 + 0.38 sp_{4\ i} = 22.56\)` - `\(\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} = 21.12\)` --- 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