class: middle, left, inverse, title-slide .title[ # Тестирование гипотез <br/>при помощи perMANOVA ] .subtitle[ ## Анализ и визуализация многомерных данных с использованием R ] .author[ ### Марина Варфоломеева ] .author[ ### Вадим Хайтов ] --- - Множественные сравнения - perMANOVA - Условия применимости perMANOVA - Более подробная интерпретация результатов perMANOVA - Более сложные дизайны в perMANOVA ### Вы сможете - Вспомнить основные дизайны для тестирования гипотез в рамках дисперсионного анализа - Применить функции R для тестирования гипотез с помощью perMANOVA --- class: middle, center, inverse # Еще раз про множественные сравнения --- ## В чем опасность множественных сравнений? -- `\(\alpha\)` --- это __вероятность совершить ошибку первого рода при тестировании гипотезы__ (= вероятность отвергнуть истинную нулевую гипотезу, = вероятность найти различия там, где их нет). Обычно принимается, что `\(H_0\)` отвергают на уровне значимости `\(\alpha = 0.05\)`. -- Когда у нас два средних --- все просто, сравнение всего одно. Естественно, __вероятность совершить ошибку I рода для группы сравнений `\(\alpha _{familywise}\)` равна уровню значимости для единственного сравнения `\(\alpha _{per\,comparison}\)`__. `\(\alpha _{familywise} = \alpha _{per\,comparison}\)` Но если сравнений много, то станет выше вероятность совершить хотя бы одну ошибку I рода (найти различия там, где их нет). --- ## Если сравнений много... Например, если мы хотим попарно сравнить три значения, нам понадобится сделать 3 сравнения. Пусть мы решили, что в каждом из сравнений будем использовать уровень значимости `\(0.05\)`. -- Тогда в каждом из сравнений вероятность совершить ошибку первого рода будет `\(\alpha_{per\,comparison} = 0.05\)`. -- .content-box-red[ Если сделать `\(N\)` тестов, то вероятность совершить хотя бы одну ошибку I рода в группе тестов (family-wise error rate, FWER) значительно возрастает. `$$\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N$$` ] --- ## Чем больше сравнений, тем больше вероятность обнаружить различия там, где их на самом деле нет. `$$\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N$$` В таблице даны значения `\(\alpha _{familywise}\)` для разного числа сравнений, если `\(\alpha _{per\,comparison} = 0.05\)`: <table> <thead> <tr> <th style="text-align:right;"> Число средних </th> <th style="text-align:right;"> Число сравнений </th> <th style="text-align:right;"> α<sub>familywise</sub> </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 0.0500 </td> </tr> <tr> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 0.1426 </td> </tr> <tr> <td style="text-align:right;"> 4 </td> <td style="text-align:right;"> 6 </td> <td style="text-align:right;"> 0.2649 </td> </tr> <tr> <td style="text-align:right;"> 5 </td> <td style="text-align:right;"> 10 </td> <td style="text-align:right;"> 0.4013 </td> </tr> </tbody> </table> --- ## Для решения проблемы есть два подхода 1. Взять более жесткий порог уровня значимости. Например, - поправка Бонферрони, - поправка Хольма-Бонферрони, - процедура Бенъямини-Хохберга. 2. Изменить схему тестирования гипотезы --- тестировать не три независимых гипотезы, а одну сложную (так это, например, происходит в ANOVA). --- ## Поправка Бонферрони `$$\alpha^*_{per\,comparison} = \alpha_{familywise} /{N}$$` .content-box-red[ Это жесткий способ, т.к с возрастанием числа сравнений резко снижается уровень значимости и мощность теста. В результате растет риск не найти различий, где они есть. ] Ниже `\(\alpha _{per\,comparison}\)` после поправки Бонферрони, сохраняющие `\(\alpha _{familywise} = 0.05\)`: <table> <thead> <tr> <th style="text-align:right;"> Число средних </th> <th style="text-align:right;"> Число сравнений </th> <th style="text-align:right;"> α<sub>per comparison</sub> </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 0.0500 </td> </tr> <tr> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 0.0167 </td> </tr> <tr> <td style="text-align:right;"> 4 </td> <td style="text-align:right;"> 6 </td> <td style="text-align:right;"> 0.0083 </td> </tr> <tr> <td style="text-align:right;"> 5 </td> <td style="text-align:right;"> 10 </td> <td style="text-align:right;"> 0.0050 </td> </tr> <tr> <td style="text-align:right;"> 6 </td> <td style="text-align:right;"> 15 </td> <td style="text-align:right;"> 0.0033 </td> </tr> </tbody> </table> --- ## Метод Хольма-Бонферрони Метод Хольма-Бонферрони --- это пошаговая процедура. Чтобы зафиксировать `\(FWER \le \alpha_{familywise}\)`: 1. Сортируем в порядке возрастания `\(N\)` значений `\(p\)`, полученные в тестах, и присваиваем каждой ранг `\(j\)` от 1 до `\(N\)`: `$$p_{1} \le p_{2} \le \cdots \le p_{N - 1} \le p_{N}$$` 2. Вводим поправку для значений `\(p\)` `$${p^*_{j}} = min{\{(N - j + 1) \cdot p_{j},\;1\}}$$` 3. Сравниваем с `\(\alpha\)`. Отвергаем `\(H_0\)`, если `\(p^*_{j} < \alpha_{familywise}\)` --- ## Пример поправки Хольма-Бонферрони В таблице ниже даны результаты нескольких сравнений (при `\(N = 5\)`). С помощью поправки Хольма-Бонферрони для каждого `\(p_j\)` мы получим свое скорректированное значение `\(p^*_{j}\)`. | Ранг(_j_)| `\(\mathbf{p_{j}}\)`| `\(N - j + 1\)`| `\(\mathbf{p^*_{j}} = min{\{(N - j + 1) \cdot p_{j},\;1\}}\)`|Отвергаем `\(H_0\)`? | |---------:|----------------:|-----------:|---------------------------------------------------------:|:-----------------------------------| | 1| 0.001| 5| 0.005|Да | | 2| 0.010| 4| 0.040|Да | | 3| 0.035| 3| 0.105|Нет, и дальше везде сохраняем `\(H_0\)` | | 4| 0.040| 2| 0.080|Нет | | 5| 0.046| 1| 0.046|Нет | --- class: middle, center, inverse # Дисперсионный анализ (повторение) --- ## Классический дисперсионный анализ .pull-left-33[  .small[Рональд Элмер Фишер] ] .pull-right-66[ Пусть имеется несколько градаций фактора __A__ (например, `\(A _{1...3}\)`) - Почему появляется __межгрупповая__ изменчивость, то есть разные средние значения для групп по фактору _А_? - Почему появляется __внутригрупповая__ изменчивость, то есть разные значения _y_ в группах? ] | A1 | A2 | A3 | | |: --------- :|: --------- :|: --------- :|: -- :| | `\(y_{11}, y_{12}, y_{13}, ...\)` | `\(y_{21}, y_{22}, y_{23}, ...\)` | `\(y_{31}, y_{32}, y_{33}, ...\)` | | | `\(\bar{y}_{A_1}\)` | `\(\bar{y}_{A_2}\)` | `\(\bar{y}_{A_3}\)` | `\(\bar{Y}_{gen}\)` | --- ## Суммарная дисперсия может быть разложена на две составляющие `\(SS_{total} = SS_{between} + SS_{within}\)` <!-- --> --- ## Тестирование гипотезы о влиянии фактора Чтобы сравнить межгрупповую изменчивость ( `\(SS_{between}\)`) и внутригрупповую ( `\(SS_{within}\)`), Фишер предложил статистику `$$F = \frac{SS_{between} / (a - 1)}{SS_{within} / (N - a)}$$` где `\(a\)` --- число групп. ### F-распределение Если межгрупповая изменчивость равна внутригрупповой ( `\(H_0\)`), то `\(F\)` принадлежит `\(F\)`-распределению с двумя параметрами `\(df_{between} = a - 1\)` и `\(df_{within} = N - a\)`, где `\(a\)` --- число классов, `\(N\)` --- общее количество объектов. <!-- --> --- ## Однофакторный дизайн Выявляется влияние фактора А.  --- ## Многофакторный ортогональный дизайн Выявляется влияние фактора А, B и их взаимодействия A `\(\times\)` B.  --- ## Иерархический дизайн Некоторые дискретные предикторы могут быть иерархически соподчинены. Выявляется влияние вложенного фактора B внутри градаций фактора A.  Влияние вложенного фактора B можно оценить, но чаще всего оно не представляет интереса для исследователя. --- ## Рандомизированный полный блочный дизайн  Влияние блокриующего фактора тоже можно оценить, но часто оно не представляет интереса для исследователя. Обычно блокирующий фактор рассматривается как случайный. --- class: middle, center, inverse # Пример: Поведение песчанок в тесте открытое поле --- class: split-20 .row.bg-main1[.content[ ## Пример: Поведение песчанок в тесте открытое поле **Гипотеза:** Разные виды песчанок различаются по поведению в тесте "Открытое поле" ]] .row.bg-main2[ .split-three[ .column[ <img src="images/Peschanka-karl.jpg" height="180" alt="Карликовая песчанка"> .small[Карликовая песчанка _Gerbillus gerbillus_] .tiny[XV8-Crisis] ] .column[ <img src="images/Peschanka-mongol.jpg" height="180" alt="Монгольская песчанка"> .small[Монгольская песчанка _Meriones unguiculatus_] .tiny[Alastair Rae] ] .column[ <img src="images/Peschanka-zhirn.jpg" height="180" alt="Жирнохвостая песчанка"> .small[Жирнохвостая песчанка _Pachyuromys duprasi_] .tiny[P.H.J. (Peter) Maas / www.petermaas.nl] ] ]] --- ## Тест "открытое поле" <div align="center"> <iframe width="560" height="500" src="https://www.youtube.com/embed/LifsadrTAUY?rel=0&start=73" frameborder="0" allowfullscreen></iframe> </div> --- class: split-20 .row.bg-main1[.content[ ## Пример: Поведение песчанок в тесте открытое поле **Гипотеза:** Разные виды песчанок различаются по поведению в тесте "Открытое поле" ]] .row.bg-main2[.content[ .split-two[ .split-three[ .column[ <img src="images/Peschanka-karl.jpg" height="180" alt="Карликовая песчанка"> .small[Карликовая песчанка _Gerbillus gerbillus_] .tiny[XV8-Crisis] ] .column[ <img src="images/Peschanka-mongol.jpg" height="180" alt="Монгольская песчанка"> .small[Монгольская песчанка _Meriones unguiculatus_] .tiny[Alastair Rae] ] .column[ <img src="images/Peschanka-zhirn.jpg" height="180" alt="Жирнохвостая песчанка"> .small[Жирнохвостая песчанка _Pachyuromys duprasi_] .tiny[P.H.J. (Peter) Maas / www.petermaas.nl] ] ] .row[ <br/> 1. Время до выхода в квадрат открытого поля 1. Количество актов мочеиспускания 1. Количество актов дефекации 1. Количество пересеченных квадратов 1. Число вертикальных стоек 1. Количество актов смещенной активности 1. Время проведенное в центре квадрата открытого поля ] ]]] --- ## Задание - Логарифмируйте исходные данные и сохраните в переменную `log_pesch`. Какие еще преобразования могли бы подойти кроме логарифма? - Какие меры различия подходят для этих данных? Почему не подходит коэффициент Брея-Куртиса? - Сделайте ординацию nMDS. Каков ее стресс? - Сохраните в датафрейм координаты точек после nMDS, добавьте в него переменную `Species` из исходных данных. - Постройте график ординации и раскрасьте точки в соответствии с видами. ```r ## Данные наблюдений pesch <- read.csv("data/pesch.csv", header = TRUE, sep = ";") head(pesch) ``` <div data-pagedtable="false"> <script data-pagedtable-source type="application/json"> {"columns":[{"label":[""],"name":["_rn_"],"type":[""],"align":["left"]},{"label":["Gender"],"name":[1],"type":["chr"],"align":["left"]},{"label":["Species"],"name":[2],"type":["chr"],"align":["left"]},{"label":["Time_to_entrance"],"name":[3],"type":["int"],"align":["right"]},{"label":["Urination"],"name":[4],"type":["int"],"align":["right"]},{"label":["Defecation"],"name":[5],"type":["int"],"align":["right"]},{"label":["Quadr_Number"],"name":[6],"type":["int"],"align":["right"]},{"label":["Vert_Number"],"name":[7],"type":["int"],"align":["right"]},{"label":["Displ_Act"],"name":[8],"type":["int"],"align":["right"]},{"label":["Time_in_Centre"],"name":[9],"type":["int"],"align":["right"]}],"data":[{"1":"f","2":"karl","3":"0","4":"1","5":"0","6":"47","7":"11","8":"4","9":"0","_rn_":"1"},{"1":"f","2":"karl","3":"20","4":"0","5":"3","6":"317","7":"58","8":"3","9":"6","_rn_":"2"},{"1":"m","2":"karl","3":"181","4":"0","5":"0","6":"177","7":"55","8":"6","9":"2","_rn_":"3"},{"1":"f","2":"karl","3":"0","4":"0","5":"0","6":"32","7":"5","8":"3","9":"0","_rn_":"4"},{"1":"f","2":"karl","3":"139","4":"0","5":"0","6":"205","7":"29","8":"10","9":"3","_rn_":"5"},{"1":"m","2":"karl","3":"0","4":"0","5":"2","6":"38","7":"9","8":"8","9":"0","_rn_":"6"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}} </script> </div> --- ## Решение Поскольку измеренные признаки варьируют в разных масштабах, целесообразно логарифмировать данные, чтобы "уменьшить" очень высокие значения. ```r options(digits = 4) log_pesch <- pesch log_pesch[, 3:ncol(pesch)] <- log(pesch[, 3:ncol(pesch)] + 1) head(log_pesch) ``` ``` Gender Species Time_to_entrance Urination Defecation Quadr_Number Vert_Number 1 f karl 0.000 0.6931 0.000 3.871 2.485 2 f karl 3.045 0.0000 1.386 5.762 4.078 3 m karl 5.204 0.0000 0.000 5.182 4.025 4 f karl 0.000 0.0000 0.000 3.497 1.792 5 f karl 4.942 0.0000 0.000 5.328 3.401 6 m karl 0.000 0.0000 1.099 3.664 2.303 Displ_Act Time_in_Centre 1 1.609 0.000 2 1.386 1.946 3 1.946 1.099 4 1.386 0.000 5 2.398 1.386 6 2.197 0.000 ``` --- ## Решение ```r library(ggplot2) library(vegan) theme_set(theme_bw(base_size = 14)) ord_pesch <- metaMDS(log_pesch[, 3:ncol(pesch)], distance = "euclidean") dfr_pesch <- as.data.frame(ord_pesch$points) dfr_pesch$Species <- pesch$Species gg_ord_pesch <- ggplot(dfr_pesch, aes(x = MDS1, y = MDS2, colour = Species)) + geom_point(size = 5) + coord_fixed() gg_ord_pesch ``` <!-- --> --- ## Различаются ли виды песчанок по поведению? Каким способом можно ответить на этот вопрос? <!-- --> -- - Мы могли бы взять каждый из поведенческих признаков отдельно и провести одномерный однофакторный дисперсионный анализ. -- - Но нас интересует поведение в целом --- нужен многомерный анализ. --- class: middle, center, inverse # Методы выявления различий между группами для многомерных данных --- ## ANOVA разработан для одномерных данных ### Что делать если мы хотим оценивать объект по многим признакам сразу? Примеры: * Сообщество как целое (много видов) * Поведение как целое (много элементов) * Метаболом как целое (много метаболитов) * Ответы респондентов на множество вопросов в анкетах -- Варианты решений: - MANOVA (Fisher, 1925, Wilks, 1932) - perMANOVA (Anderson, 2001; McArdle, Anderson, 2001) - distance-based Redundancy Analysis (db-RDA) (Legendre, Anderson, 1999) --- ## Многомерный дисперсионный анализ <br/>(MANOVA, Multivariate Analysis Of Variance) - Разработан давно - Аналог ANOVA для многомерных данных - Многомерное нормальное распределение данных В MANOVA сравниваются: - __.red[отклонения точек от групповых центроидов]__ (аналог `\(SS_{within}\)`) - __.blue[отклонения групповых центроидов от общего центроида]__ (аналог `\(SS_{between}\)`). .pull-left[ <img src="images/MANOVA_mod.png" height="260" alt="MANOVA"> .small[Anderson, 2001] ] -- .pull-right[ Ограничения MANOVA: - Многомерная нормальность распределения - Гомогенность дисперсий ] --- ## Permutational Multivariate Analysis of Variance (perMANOVA) .pull-left-33[  .small[Марти Джейн Андерсон] ] .pull-right-66[  ] --- ## Теорема .content-box-grey[ Сумма квадратов __.red[отклонений от центроидов]__ равна сумме квадратов __.orange[взаимных расстояний]__, деленной на число объектов ] Для Евклидовых расстояний эта закономерность была известна давно (например, Kendall, Stuart 1963). .pull-left[  .tiny[Anderson, 2001] ] -- .pull-right[ В случае Евклидова расстояния (именно его имплицитно использует MANOVA) центроиды найти просто --- это средние соответствующих координат. Поэтому обычно сначала вычисляли центроиды, а потом --- сумму квадратов отклонений от них. ] -- Для других мер различия центроиды найти гораздо сложнее. Например, для коэффициента Брея-Куртиса (не метрика), среднее значение координат не будет соответствовать центроиду. --- ## Марти Андерсон показала, что можно обойтись без вычисления центроидов и для других мер различия <center><img src="images/Centroid and distance.png" height="170" alt="Centroid and distance"></center><small>Anderson, 2001</small> .pull-left[ <br/> MANOVA (Евклидово расстояние) <img src="images/MANOVA.png" height="260" alt="MANOVA"><small>Anderson, 2001</small> ] .pull-right[ perMANOVA (любой коэффициент) <img src="images/distances-permanova.png" height="260" alt="perMANOVA"> ] --- ## Можно непосредственно из матрицы любых коэффициентов различия найти и общую и внутригрупповые суммы квадратов <img src="images/Matrix transformation.png" height="500" alt=" "><small>Anderson, 2001</small> --- ## Разложение дисперсии становится очень простым .pull-left[ <img src="images/distances-total-within-between.png" width="450" alt="distances-total-within-between"> <small>Anderson, 2001</small> ] .pull-right[ Пусть всего _N_ элементов, принадлежащих _a_ группам по _n_ элементов в каждой, _d_ - расстояние между _i_-тым и _j_-тым объектами, `\(\epsilon\)` - 1 если объекты _i_ и _j_ из одной группы и 0, если из разных. ] `\(SS_{total} = \frac{1}{N}\sum \sum {d_{ij}^2}\)` Сумма квадратов взаимных расстояний --- это сумма квадратов субдиагональных элементов, деленная на число объектов _N_. `\(SS_{within} = \frac{1}{n}\sum \sum {d_{ij}^2 \cdot \epsilon_{ij}}\)` Внутригрупповая сумма квадратов --- это сумма всех сумм квадратов расстояний между элементами для каждой группы, деленная на _n_ число объектов в группе. Тогда межгрупповая сумма квадратов `\(SS_{between} = SS_{total} - SS_{within}\)` --- ## Псевдо-F статистика `$$F = \frac{SS_{between} / (a - 1)}{SS_{within}/(N - a)}$$` Для оценки значимости псевдо-F используется __пермутационная процедура__: - Случайным образом перетасовываются строки исходной матрицы - После каждого акта пермутации вычисляется `\(F_{perm}\)` - Уровень значимости `$$p = \frac {N_{F_{perm} \ge F}} {N_{permutations}}$$` -- .content-box-red[ Внимание! Для иерархического дизайна процедура пермутации имеет свои особенности (обсудим позднее). ] --- class: middle, center, inverse # perMANOVA в примере --- ## Применим метод perMANOVA (функция `adonis()`) ```r library(vegan) permanova_pesch <- adonis2(log_pesch[3:9] ~ Species, data = log_pesch, method = "euclidean") permanova_pesch ``` ``` Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 adonis2(formula = log_pesch[3:9] ~ Species, data = log_pesch, method = "euclidean") Df SumOfSqs R2 F Pr(>F) Species 2 60 0.237 3.58 0.012 * Residual 23 193 0.763 Total 25 253 1.000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Мы видим традиционную для ANOVA таблицу результатов. Что здесь что? --- ## Результаты perMANOVA ``` Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 adonis2(formula = log_pesch[3:9] ~ Species, data = log_pesch, method = "euclidean") Df SumOfSqs R2 F Pr(>F) Species 2 60 0.237 3.58 0.012 * Residual 23 193 0.763 Total 25 253 1.000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` -- - Поведение разных видов песчанок в тесте открытое поле значимо различалось (perMANOVA, p = 0.012). --- class: middle, center, inverse # Условия применимости perMANOVA --- ## Условия применимости perMANOVA .content-box-purple[ 1. Требуется равенство внутригрупповых дисперсий (гомоскедастичность). 2. Желательно использование сбалансированных данных (с равным объемом групп), т.к. этом случае perMANOVA устойчив к гетерогенности дисперсии (Anderson, Walsh, 2013). ] --- ## Проверка равенства внутригрупповых дисперсий Для проверки можно использовать функцию `betadisper()`, изначально предназначенную для сравнения `\(\beta\)`-разнообразия сообществ. -- Функция `betadisper()` вычисляет внутригрупповые центроиды и координаты точек в пространстве главных координат (_Principal Coordinate Analysis = PCoA = metric MDS_). Затем, значимость различий отклонений от центроидов в разных группах проверяется с помощью процедуры `PERMDISP2`. <div id="htmlwidget-5d15f6ffc3281a939cbc" style="width:850px;height:180px;" class="grViz html-widget"></div> <script type="application/json" data-for="htmlwidget-5d15f6ffc3281a939cbc">{"x":{"diagram":"digraph {\n graph [layout = dot, rankdir = LR, ratio = compress]\n node [shape = rectangle, fontsize = 16, fontname=\"Helvetica\", fixedsize = true, height = 1.5, style = filled, color=\"#E8E8E8\"];\n edge [ fontname=Courier, fontsize = 18];\n rec1 [label = \"Исходные \\nданные\", width = 1.5]\n rec2 [label = \"Матрица \\nкоэффициентов различия, \\nв т.ч. неевклидовых\", width = 3]\n rec3 [label = \"PCoA: \\nцентроиды\\nи координаты точек, \\nЕвклидово пространство\", width = 3]\n# rec4 [label = \"График \\nординации\\nPCoA\", width = 2]\n rec5 [label = \"PERMDISP2\\n Тест различий \\nвнутригрупповых \\nдисперсий\", width = 2]\n # edge definitions with the node IDs\n rec1 -> rec2 [ label=\"vegdist()\" ];\n rec2 -> rec3 [label = \"betadisper()\"];\n# rec3 -> rec4 [label = \"plot()\"];\n rec3 -> rec5 [label = \"anova()\"];\n\n }","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script> <br/> -- ```r dist_pesch <- vegdist(log_pesch[,3:ncol(pesch)], method = "euclidean") PCO_pesch <- betadisper(dist_pesch, log_pesch$Species) ``` --- ## График ординации PCoA Объект, возвращаемый `betadisper()`, позволяет также нарисовать наши объекты в пространстве главных координат (PCoA). ```r plot(PCO_pesch, main = "PCoA ordination") ``` <!-- --> --- ## Процедура `PERMDISP2` для проверки равенства внутригрупповых дисперсий Процедура `PERMDISP2` реализована в пакете `vegan` в функции `anova.betadisper()`. Это многомерный аналог теста Левина на гомогенность дисперсий в группах, который иногда используется для проверки условий применимости дисперсионного анализа. -- ```r anova(PCO_pesch) ``` ``` Analysis of Variance Table Response: Distances Df Sum Sq Mean Sq F value Pr(>F) Groups 2 9.7 4.86 2.05 0.15 Residuals 23 54.5 2.37 ``` - Не выявлено значимых различий разброса внутригрупповых расстояний. --- ## Для визуализации можно нарисовать боксплот ```r boxplot(PCO_pesch) ``` <!-- --> --- class: middle, center, inverse # Более подробная интерпретация результатов perMANOVA --- ## Post hoc тесты в perMANOVA <!-- --> На приведенной ординации видно, что точки, соответствующие Монгольским песчанкам расположены отдельно от остальных. Для выявления попарных различий нужны попарные сравнения. --- class: split-20 .row.bg-main1[.content[ ## Попарные сравнения как замена post hoc тесту Внимание! В пакете `vegan` пост хок тест не реализован. Но мы можем сделать простейшую его версию самостоятельно. ]] .row.bg-main2[.content[ .split-two[ .split-three[ .column[ <img src="images/Peschanka-karl.jpg" height="180" alt="Карликовая песчанка"> .small[Карликовая песчанка _Gerbillus gerbillus_] .tiny[XV8-Crisis] ] .column[ <img src="images/Peschanka-mongol.jpg" height="180" alt="Монгольская песчанка"> .small[Монгольская песчанка _Meriones unguiculatus_] .tiny[Alastair Rae] ] .column[ <img src="images/Peschanka-zhirn.jpg" height="180" alt="Жирнохвостая песчанка"> .small[Жирнохвостая песчанка _Pachyuromys duprasi_] .tiny[P.H.J. (Peter) Maas / www.petermaas.nl] ] ] .row[ <br/> Проведем __попарные сравнения__ между группами, то есть: * Карликовые _vs._ Монгольские * Карликовые _vs._ Жирнохвостые * Монгольские _vs._ Жирнохвостые ] ]]] --- ## Функция для попарных сравнений perMANOVA ```r pairwise_permanova <- function(dat, group, strata = NULL, ...){ groups <- unique(as.character(group)) all_pairs <- combn(groups, 2, simplify = FALSE) ncomb <- length(all_pairs) res <- rep(NA, ncomb) for (i in 1:ncomb) { flt <- group %in% all_pairs[[i]] strata_flt <- if(is.null(strata)) strata[flt] posthoc <- adonis2(dat[flt, ] ~ group[flt], strata = strata_flt, ...) res[i] <- posthoc$`Pr(>F)`[1] names(res)[i] <- paste0(all_pairs[[i]][1], "_", all_pairs[[i]][2]) } return(res) } ``` --- ## Результаты попарных сравнений .center[    ] ```r p_vals <- pairwise_permanova( dat = log_pesch[, -c(1:2)], group = log_pesch$Species, method = "euclidean", permutations=9999) p_vals ``` ``` karl_mongol karl_zhirnokhvost mongol_zhirnokhvost 0.0011 0.4029 0.0031 ``` Это все? Пишем статью? <br/><br/> .small[Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl ] --- ## Поправка на множественные сравнения Мы делали три пары сравнений --- нужно ввести поправку для уровня значимости. Будем считать значимыми различия в тех сравнениях, где после введения поправки `\(p < 0.05\)`. Можно посчитать скорректированные значения уровня значимости `\(p\)` с учетом поправки Хольма-Бонферрони. ```r p_vals # было без поправки ``` ``` karl_mongol karl_zhirnokhvost mongol_zhirnokhvost 0.0011 0.4029 0.0031 ``` ```r p.adjust(p_vals, method = "holm") # стало с поправкой ``` ``` karl_mongol karl_zhirnokhvost mongol_zhirnokhvost 0.0033 0.4029 0.0062 ``` --- class: middle, center, inverse # Более сложные дизайны в perMANOVA --- ## Многофакторный дизайн в perMANOVA Выясним, влияет ли пол и вид песчанок на поведение. Отфильтруем исходные данные (в случае с жирнохвостыми песчанками были изучены только самки) ```r log_pesch2 <- log_pesch[log_pesch$Species != "zhirnokhvost", ] ``` --- ## Проведем двухфакторный анализ perMANOVA Различается ли поведение песчанок в зависимости от видовой принадлежности и пола? ```r twofact_pesch <- adonis2(log_pesch2[,3:ncol(pesch)] ~ Gender * Species, data = log_pesch2, method = "euclidean") twofact_pesch ``` ``` Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 adonis2(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender * Species, data = log_pesch2, method = "euclidean") Df SumOfSqs R2 F Pr(>F) Gender 1 7.2 0.049 1.13 0.292 Species 1 45.5 0.311 7.16 0.004 ** Gender:Species 1 4.5 0.031 0.72 0.529 Residual 14 88.9 0.608 Total 17 146.1 1.000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Здесь возможен иерархический дизайн Различается ли поведение самцов и самок у этих видов песчанок? ```r nested_pesch <- adonis2(log_pesch2[, 3:ncol(pesch)] ~ Gender, data = log_pesch2, strata = log_pesch2$Species, method = "euclidean") nested_pesch ``` ``` Permutation test for adonis under reduced model Terms added sequentially (first to last) Blocks: strata Permutation: free Number of permutations: 999 adonis2(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender, data = log_pesch2, method = "euclidean", strata = log_pesch2$Species) Df SumOfSqs R2 F Pr(>F) Gender 1 7.2 0.049 0.83 0.3 Residual 16 138.9 0.951 Total 17 146.1 1.000 ``` Внимание! Пермутации производятся только в пределах группирующего фактора (аргумент `strata`) --- ## Задание + Создайте датафрейм из файла `simulated_data.csv` (Это данные симулированные по алгоритму, приведенному в справке по функции `adonis()`) + В этом датафрейме записано обилие двух видов на экспериментальных площадках двух типов: без добавления и с добавлением NO3, по 6 повторностей в каждом эксперименте. Эксперименты были независимо проведены на 3 полях. + Оцените, зависит ли структура совместного поселения этих двух видов от концентрации NO3. --- ## Решение ```r com <- read.csv("data/simulated_data.csv", sep = ',', header = T) adonis2(com[,1:2] ~ com$NO3) # Ошибочный дизайн ``` ``` Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 999 adonis2(formula = com[, 1:2] ~ com$NO3) Df SumOfSqs R2 F Pr(>F) com$NO3 1 0.048 0.084 3.13 0.053 . Residual 34 0.521 0.916 Total 35 0.569 1.000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r adonis2(com[,1:2] ~ com$NO3, strata = com$field) # Правильный дизайн ``` ``` Permutation test for adonis under reduced model Terms added sequentially (first to last) Blocks: strata Permutation: free Number of permutations: 999 adonis2(formula = com[, 1:2] ~ com$NO3, strata = com$field) Df SumOfSqs R2 F Pr(>F) com$NO3 1 0.048 0.084 3.13 0.001 *** Residual 34 0.521 0.916 Total 35 0.569 1.000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Summary + При одновременном тестировании нескольких гипотез растет вероятность ошибки I рода. + Чтобы контролировать вероятность ошибки I рода, либо используют более жесткий уровень значимости, либо тестируют сложную гипотезу вместо нескольких простых. + perMANOVA дает возможность тестировать сложные гипотезы в отношении явлений, описанных по многим переменным (т.е. на многомерных данных). + В perMANOVA можно использовать любые коэффициенты различия. + Для применения perMANOVA требуется равенство разбросов точек между центроидами их групп, но при равных объемах групп анализ устойчив к отклонениям от этого условия. + При использовании perMANOVA важно не запутаться в дизайне. --- ## Другие программы * *Primer 6.0 + PERMANOVA* Коммерческий продукт. * *PAST* Некоммерческая программа. Здесь метод называется NPMANOVA. * Оригинальная программа М. Андерсон *PERMANOVA* и *PERMDISP*. --- ## Что почитать * Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26: 32–46. * Anderson, M.J. 2005. PERMANOVA: a FORTRAN computer program for permutational multivariate analysis of variance. Department of Statistics, University of Auckland, New Zealand. * Anderson, M.J. (2004). PERMDISP: a FORTRAN computer program for permutational analysis of multivariate dispersions (for any two-factor ANOVA design) using permutation tests. Department of Statistics, University of Auckland, New Zealand. * Legendre P., Legendre L. (2012) Numerical ecology. Second english edition. Elsevier, Amsterdam.