Тестирование гипотез при помощи perMANOVA

# Тестирование гипотез <br/>при помощи perMANOVA
## Анализ и визуализация многомерных данных с использованием R
### Марина Варфоломеева
### Вадим Хайтов

---

- Еще раз про множественные сравнения
- Пример: Поведение песчанок в тесте открытое поле
- perMANOVA
- Условия применимости perMANOVA
- Более подробная интерпретация результатов perMANOVA
- Более сложные дизайны в perMANOVA

### Вы сможете

- Вспомнить основные дизайны для тестирования гипотез в рамках дисперсионного анализа
- Применить функции R для тестирования гипотез с помощью perMANOVA

---

# Еще раз про множественные сравнения

---

## В чем опасность множественных сравнений?

`$\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$`.

Если сделать `$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$`:

| Число средних| Число сравнений| &alpha;<sub>familywise</sub>|
|-------------:|---------------:|----------------------------:|
|             2|               1|                        0.050|
|             3|               3|                        0.143|
|             4|               6|                        0.265|
|             5|              10|                        0.401|

---

## Для решения проблемы есть два подхода

1. Взять более жесткий порог уровня значимости. Например,
  - поправка Бонферрони,
  - поправка Хольма-Бонферрони,
  - процедура Бенъямини-Хохберга.
2. Изменить схему тестирования гипотезы --- тестировать не три независимых гипотезы, а одну сложную (так это, например, происходит в ANOVA).

---

## Поправка Бонферрони

`$\alpha^*_{per\,comparison} = \alpha_{familywise} /{N}$`

Это жесткий способ, т.к с возрастанием числа сравнений резко снижается уровень значимости и мощность теста. В результате растет риск не найти различий, где они есть.

Ниже `$\alpha _{per\,comparison}$` после поправки Бонферрони, сохраняющие `$\alpha _{familywise} = 0.05$`:

| Число средних| Число сравнений| &alpha;<sub>per comparison</sub>|
|-------------:|---------------:|--------------------------------:|
|             2|               1|                            0.050|
|             3|               3|                            0.017|
|             4|               6|                            0.008|
|             5|              10|                            0.005|
|             6|              15|                            0.003|

---

## Метод Хольма-Бонферрони

Метод Хольма-Бонферрони --- это пошаговая процедура.

Чтобы зафиксировать `$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|Нет                                 |

---

# Дисперсионный анализ (повторение)

---

## Классический дисперсионный анализ

<small>Рональд Элмер Фишер</small>
]

Пусть имеется несколько градаций фактора _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$` --- общее количество объектов в анализе.

---

## Однофакторный дизайн

Выявляется влияние фактора А.

![](images/anova-designs-one-way.png)
---

## Многофакторный ортогональный дизайн

Выявляется влияние фактора А, B и их взаимодействия A `$\times$` B.

![](images/anova-designs-two-way.png)
---

## Иерархический дизайн

Некоторые дискретные предикторы могут быть иерархически соподчинены.

Выявляется влияние вложенного фактора B внутри градаций фактора A.

![](images/anova-designs-nested.png)

Влияние вложенного фактора B можно оценить, но чаще всего оно не представляет интереса для исследователя.

---

## Рандомизированный полный блочный дизайн

![](images/anova-designs-block.png)

Влияние блокриующего фактора тоже можно оценить, но часто оно не представляет интереса для исследователя. Обычно блокирующий фактор рассматривается как случайный.

---

# Пример: Поведение песчанок в тесте открытое поле

---

.row.bg-main1[.content[

## Пример: Поведение песчанок в тесте открытое поле

**Гипотеза:** Разные виды песчанок различаются по поведению в тесте "Открытое поле"

]]

.row.bg-main2[
.split-three[
.column[
<img src="images/Peschanka-karl.jpg" height="180" alt="Карликовая песчанка">

Карликовая песчанка (_Gerbillus gerbillus_)

.small[XV8-Crisis]
]
.column[
<img src="images/Peschanka-mongol.jpg" height="180" alt="Монгольская песчанка">

Монгольская песчанка (_Meriones unguiculatus_)

.small[Alastair Rae]
]
.column[
<img src="images/Peschanka-zhirn.jpg" height="180" alt="Жирнохвостая песчанка">

Жирнохвостая песчанка (_Pachyuromys duprasi_)

---

## Тест "открытое поле"

---

.row.bg-main1[.content[

## Пример: Поведение песчанок в тесте открытое поле

**Гипотеза:** Разные виды песчанок различаются по поведению в тесте "Открытое поле"

]]

.row.bg-main2[.content[
.split-two[
.split-three[
.column[
<img src="images/Peschanka-karl.jpg" height="180" alt="Карликовая песчанка">

Карликовая песчанка (_Gerbillus gerbillus_)

.small[XV8-Crisis]
]
.column[
<img src="images/Peschanka-mongol.jpg" height="180" alt="Монгольская песчанка">

Монгольская песчанка (_Meriones unguiculatus_)

.small[Alastair Rae]
]
.column[
<img src="images/Peschanka-zhirn.jpg" height="180" alt="Жирнохвостая песчанка">

Жирнохвостая песчанка (_Pachyuromys duprasi_)

.row[
<br/>
1. Время до выхода в квадрат открытого поля
1. Количество актов мочеиспускания
1. Количество актов дефекации
1. Количество пересеченных квадратов
1. Число вертикальных стоек
1. Количество актов смещенной активности
1. Время проведенное в центре квадрата открытого поля
]
]]]

---

## Задание

- Какую меру различия можно использовать с этими данными?
- Как можно преобразовать данные?
- Постройте ординацию объектов в осях MDS и раскрасьте точки в соответствии с видами

```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
## 1      f    karl            0.000    0.6931      0.000        3.871
## 2      f    karl            3.045    0.0000      1.386        5.762
## 3      m    karl            5.204    0.0000      0.000        5.182
## 4      f    karl            0.000    0.0000      0.000        3.497
## 5      f    karl            4.942    0.0000      0.000        5.328
## 6      m    karl            0.000    0.0000      1.099        3.664
##   Vert_Number Displ_Act Time_in_Centre
## 1       2.485     1.609          0.000
## 2       4.078     1.386          1.946
## 3       4.025     1.946          1.099
## 4       1.792     1.386          0.000
## 5       3.401     2.398          1.386
## 6       2.303     2.197          0.000
```

---

## Решение

```r
library(ggplot2)
library(vegan)
theme_set(theme_bw(base_size = 14))

mds_pesch <- metaMDS(log_pesch[, 3:ncol(pesch)], distance = "euclidean")
mds_pesch <- as.data.frame(mds_pesch$points)
mds_pesch$Species <- pesch$Species

ggplot(mds_pesch, aes(x = MDS1, y = MDS2, colour = Species)) + 
  geom_point(size = 5)
```

---

## Различаются ли виды песчанок по поведению?

Каким способом можно ответить на этот вопрос?

- Мы могли бы взять каждый из поведенческих признаков отдельно и провести одномерный однофакторный дисперсионный анализ.

- Но нас интересует поведение в целом --- нужен многомерный анализ.

---

# Методы выявления различий между группами для многомерных данных

---

## ANOVA разработан для одномерных данных

### Что делать если мы хотим оценивать объект по многим признакам сразу?

Примеры:

* Сообщество как целое (много видов)
* Поведение как целое (много элементов)
* Метаболом как целое (много метаболитов)
* Ответы респондентов на множество вопросов в анкетах

Варианты решений:

- MANOVA (Fisher, 1925, Wilks, 1932)
- perMANOVA (Anderson, 2001; McArdle, Anderson, 2001)
- distance-based Redundancy Analysis (db-RDA) (Legendre, Anderson, 1999)

---

## Многомерный дисперсионный анализ (MANOVA)

Давно разработан параметрический метод MANOVA (Multivariate Analysis Of Variance). Он дает возможность проводить анализ аналогичный ANOVA. В основе MANOVA лежат представление о многомерном нормальном распределении и расстояниях между центроидами.

В MANOVA сравниваются:

- __.red[отклонения точек от групповых центроидов]__ (аналог `$SS_{within}$`)
- __.blue[отклонения групповых центроидов от общего центроида]__ (аналог `$SS_{between}$`).

]

Ограничения MANOVA:

- Многомерная нормальность распределения
- Гомогенность дисперсий

]

---

## Permutational Multivariate Analysis of Variance (perMANOVA)

<small>Марти Джейн Андерсон</small>

]
.pull-right-66[
<img src="images/Anderson 2001.jpg" width="500" alt="Anderson 2001">
]

---

## Теорема

Сумма квадратов <span class="red">отклонений от центроидов</span> равна сумме квадратов <span class="orange">взаимных расстояний</span>, деленной на число объектов

Для Евклидовых расстояний эта закономерность была известна давно (например, Kendall, Stuart 1963).

<small>Anderson, 2001</small>

]

В случае Евклидова расстояния (именно его имплицитно использует MANOVA) центроиды найти очень просто --- это средние значения соответствующих координат. Поэтому обычно сначала непосредственно вычисляли центроиды, и затем --- сумму квадратов отклонений от них.

]

Для других мер различия центроиды найти гораздо сложнее. Например, для коэффициента Брея-Куртиса (не метрика), среднее значение не будет соответствовать центроиду.

---

## Марти Андерсон показала, что можно обойтись без вычисления центроидов и для других мер различия

<center><img src="images/Centroid and distance.png" height="170" alt="Centroid and distance"></center><small>Anderson, 2001</small>

<br/>
MANOVA (Евклидово расстояние)

<img src="images/MANOVA.png" height="260" alt="MANOVA"><small>Anderson, 2001</small>

]

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}}$$`
---

# perMANOVA в примере

---

## Применим метод perMANOVA (функция `adonis()`)

```r
library(vegan)
permanova_pesch <- adonis(log_pesch[3:9] ~ Species, data = log_pesch,
                          method = "euclidean")
permanova_pesch
```

```
## 
## Call:
## adonis(formula = log_pesch[3:9] ~ Species, data = log_pesch,      method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)  
## Species    2        60   30.02    3.58 0.237  0.012 *
## Residuals 23       193    8.39         0.763         
## Total     25       253                 1.000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Мы видим традиционную для ANOVA таблицу результатов. Что здесь что?

---

## Результаты perMANOVA

- Поведение разных видов песчанок в тесте открытое поле значимо различалось (perMANOVA, p = 0.012).

---

# Условия применимости perMANOVA

---

## Условия применимости perMANOVA

1. Требуется равенство внутригрупповых дисперсий (гомоскедастичность).

2. Желательно использование сбалансированных данных (с равным объемом групп), т.к. этом случае perMANOVA устойчив к гетерогенности дисперсии (Anderson, Walsh, 2013).

---

## Проверка равенства внутригрупповых дисперсий

Для проверки можно использовать функцию `betadisper()`, изначально предназначенную для сравнения `$\beta$`-разнообразия сообществ.

Функция `betadisper()` вычисляет внутригрупповые центроиды и координаты точек в пространстве главных координат (_Principal Coordinate Analysis = PCoA = metric MDS_).

Затем, значимость различий отклонений от центроидов в разных группах проверяется с помощью процедуры `PERMDISP2`.

<div id="htmlwidget-5d15f6ffc3281a939cbc" style="width:850px;height:144px;" 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)
```

<img src="04_perMANOVA_files/figure-html/unnamed-chunk-13-1.png" style="display: block; margin: auto;" />
---

# Более подробная интерпретация результатов perMANOVA

---

## Post hoc тесты в perMANOVA

На приведенной ординации видно, что точки, соответствующие Монгольским песчанкам расположены отдельно от остальных.

Для выявления попарных различий нужны попарные сравнения.

---

.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="Карликовая песчанка">

Карликовая песчанка (_Gerbillus gerbillus_)

.small[XV8-Crisis]
]
.column[
<img src="images/Peschanka-mongol.jpg" height="180" alt="Монгольская песчанка">

Монгольская песчанка (_Meriones unguiculatus_)

.small[Alastair Rae]
]
.column[
<img src="images/Peschanka-zhirn.jpg" height="180" alt="Жирнохвостая песчанка">

Жирнохвостая песчанка (_Pachyuromys duprasi_)

* Карликовые _vs._ Монгольские
* Карликовые _vs._ Жирнохвостые
* Монгольские _vs._ Жирнохвостые
]
]]]

---

## Функция для попарных сравнений perMANOVA

```r
pairwise_permanova <- function(dat, group, strata = NULL, ...){
  pair <- combn(unique(as.character(group)), 2)
  ncomb <- ncol(pair)
  res <- rep(NA, ncomb)
  for (i in 1:ncomb) {
    filter <- group %in% pair[, i]
    if(is.null(strata)){
      posthoc <- adonis(dat[filter, ] ~ group[filter], ...)$aov.tab$Pr[1]
    } else {
      posthoc <- adonis(dat[filter, ] ~ group[filter], 
                        strata = strata[filter], ...)$aov.tab$Pr[1]
    }
    res[i] <- posthoc
    names(res)[i] <- paste(pair[, i], collapse = " vs. ")
  }
  return(res)
}
```

---

## Результаты попарных сравнений

```r
p_vals <- pairwise_permanova(
  dat = log_pesch[, -c(1:2)], group = log_pesch$Species, 
  method = "euclidean", permutations=9999)
p_vals
```

```
##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. 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 </small>

---

## Поправка на множественные сравнения

Мы делали три пары сравнений --- нужно ввести поправку для уровня значимости. Будем считать значимыми различия в тех сравнениях, где после введения поправки `$p < 0.05$`.

Можно посчитать скорректированные значения уровня значимости `$p$` с учетом поправки Хольма-Бонферрони.

```r
p.adjust(p_vals, method = "holm")
```

```
##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. zhirnokhvost 
##                  0.0033                  0.4029                  0.0062
```

---

# Более сложные дизайны в perMANOVA

---

## Многофакторный дизайн в perMANOVA
 
Выясним, влияет ли пол и вид песчанок на поведение.

Отфильтруем исходные данные (в случае с жирнохвостыми песчанками были изучены только самки)

```r
log_pesch2 <- log_pesch[log_pesch$Species != "zhirnokhvost", ]
```

---

## Проведем двухфакторный анализ perMANOVA

Различается ли поведение песчанок в зависимости от видовой принадлежности и пола?

```r
twofact_pesch <- adonis(log_pesch2[,3:ncol(pesch)] ~ Gender * Species,
                        data = log_pesch2, method = "euclidean")
twofact_pesch
```

```
## 
## Call:
## adonis(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender * Species,      data = log_pesch2, method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)   
## Gender          1       7.2     7.2    1.13 0.049  0.292   
## Species         1      45.5    45.5    7.16 0.311  0.004 **
## Gender:Species  1       4.5     4.5    0.72 0.031  0.529   
## Residuals      14      88.9     6.4         0.608          
## Total          17     146.1                 1.000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Здесь возможен иерархический дизайн

Различается ли поведение самцов и самок у этих видов песчанок?

```r
nested_pesch <- adonis(log_pesch2[, 3:ncol(pesch)] ~ Gender, 
                       data = log_pesch2, strata = log_pesch2$Species, 
                       method = "euclidean")
nested_pesch
```

```
## 
## Call:
## adonis(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender, data = log_pesch2,      method = "euclidean", strata = log_pesch2$Species) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)
## Gender     1       7.2    7.21    0.83 0.049    0.3
## Residuals 16     138.9    8.68         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)

# Ошибочный дизай
com_permanova <- adonis(com[,1:2] ~ com$NO3)

# Правильный дизайн
com_permanova2 <- adonis(com[,1:2] ~ com$NO3, strata = com$field)
```

---

## 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.