Неметрическое многомерное шкалирование, envfit, ordisurf

.title[
# Неметрическое многомерное шкалирование, envfit, ordisurf
]
.subtitle[
## Анализ и визуализация многомерных данных с использованием R
]
.author[
### Марина Варфоломеева
]
.author[
### Вадим Хайтов
]
.author[
### Анастасия Лянгузова
]

---

- Неметрическое многомерное шкалирование
- Как работает неметрическое многомерное шкалирование
- Оценка качества подгонки ординации
- Трактовка результатов ординации

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

- Построить диаграмму nMDS.
- Охарактеризовать качество ординации с помощью величины стресса.
- Сравнить результаты нескольких ординаций.

---

# Проблема изображения многомерных данных

---

## Облако точек в многомерном пространстве

Когда признаков много, можно представить все объекты как облако точек в многомерном пространстве.

![](./images/Migration-DonMcCullough-Flickr.jpg)

Migration by Don McCullough
on [Flickr](https://flic.kr/p/fEFhCj)

---

## Для изображения N объектов в идеале нужно N-1 измерений

- 2 объекта = 1D прямая
- 3 объекта = 2D плоскость
- 4 объекта = 3D объем
- ...
- N объектов = (N-1)-мерное пространство

---

## Многомерное пространство сложно изобразить. <br/> Есть два пути:

- Выбрать наиболее информативную проекцию (не все видно).
- Сохранить отношения между объектами (придется исказить расстояния).

.small[
black shadows for a white horses / les negres ombres dels cavalls blancs by  Ferran Jordà
on [Flickr](https://flic.kr/p/9XJxiL)
]

---

## Многомерное шкалирование (MDS) и неметрическое многомерное шкалирование (nMDS)

- Используют разложение матриц на собственные числа и вектора (про матричные операции --- чуть дальше в курсе).
]

.pull-right[
- Оси в полученной ординации случайны, и изменчивость не максимизируется вдоль осей. 
- Про алгоритм nMDS и будем сегодня говорить :) 
]

Исходные данные --- матрица расстояний между объектами в многомерном пространстве.

---

## Неметрическое многомерное шкалирование

__Nonmetric Multidimensional Scaling__, __nMDS__ (Shepard 1962, 1966, Kruskal 1964a, b)

Метод визуализации отношений между объектами в пространстве с небольшим числом измерений. Число измерений задаётся исследователем (обычно 2, реже 3).

nMDS ординация старается сохранить отношения между объектами: близкие объекты будут сближены, отдалённые отдалены друг от друга.

---

## Неметрическое многомерное шкалирование

Таким образом, зная расстояния между городами Европы,

```
          Athens Barcelona Brussels Calais
Barcelona   3313                          
Brussels    2963      1318                
Calais      3175      1326      204       
Cherbourg   3339      1294      583    460
```

мы бы смогли восстановить по ним карту.

![](02_nMDS_files/figure-html/eu-1.png)

---

## Карта пригородов Санкт-Петербурга

![](02_nMDS_files/figure-html/gg-spb-1.png)

---

## Расстояния по автодорогам

<div style="border: 1px solid #ddd; padding: 0px; overflow-y: scroll; height:500px; overflow-x: scroll; width:100%; "><table class="table" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;position: sticky; top:0; background-color: #FFFFFF;"> name </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Всеволожск </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Гатчина </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Зеленогорск </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Красное Село </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Кронштадт </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Павловск </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Парголово </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Петергоф </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Пушкин </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Санкт-Петербург </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Сертолово </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Сестрорецк </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Стрельна </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> Шлиссельбург </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 71 </td>
   <td style="text-align:right;"> 71 </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 57 </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 73 </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 24 </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 56 </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 40 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:right;"> 71 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 97 </td>
   <td style="text-align:right;"> 20 </td>
   <td style="text-align:right;"> 76 </td>
   <td style="text-align:right;"> 30 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 47 </td>
   <td style="text-align:right;"> 29 </td>
   <td style="text-align:right;"> 47 </td>
   <td style="text-align:right;"> 79 </td>
   <td style="text-align:right;"> 82 </td>
   <td style="text-align:right;"> 38 </td>
   <td style="text-align:right;"> 80 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:right;"> 71 </td>
   <td style="text-align:right;"> 97 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 81 </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 87 </td>
   <td style="text-align:right;"> 44 </td>
   <td style="text-align:right;"> 77 </td>
   <td style="text-align:right;"> 82 </td>
   <td style="text-align:right;"> 53 </td>
   <td style="text-align:right;"> 33 </td>
   <td style="text-align:right;"> 25 </td>
   <td style="text-align:right;"> 80 </td>
   <td style="text-align:right;"> 106 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 20 </td>
   <td style="text-align:right;"> 81 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 56 </td>
   <td style="text-align:right;"> 24 </td>
   <td style="text-align:right;"> 53 </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 21 </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 62 </td>
   <td style="text-align:right;"> 65 </td>
   <td style="text-align:right;"> 19 </td>
   <td style="text-align:right;"> 70 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 76 </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 56 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 75 </td>
   <td style="text-align:right;"> 42 </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 51 </td>
   <td style="text-align:right;"> 45 </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 105 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:right;"> 57 </td>
   <td style="text-align:right;"> 30 </td>
   <td style="text-align:right;"> 87 </td>
   <td style="text-align:right;"> 24 </td>
   <td style="text-align:right;"> 75 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 63 </td>
   <td style="text-align:right;"> 43 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 35 </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 34 </td>
   <td style="text-align:right;"> 64 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 44 </td>
   <td style="text-align:right;"> 53 </td>
   <td style="text-align:right;"> 42 </td>
   <td style="text-align:right;"> 63 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 68 </td>
   <td style="text-align:right;"> 59 </td>
   <td style="text-align:right;"> 18 </td>
   <td style="text-align:right;"> 9 </td>
   <td style="text-align:right;"> 29 </td>
   <td style="text-align:right;"> 54 </td>
   <td style="text-align:right;"> 66 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 73 </td>
   <td style="text-align:right;"> 47 </td>
   <td style="text-align:right;"> 77 </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 43 </td>
   <td style="text-align:right;"> 68 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 41 </td>
   <td style="text-align:right;"> 72 </td>
   <td style="text-align:right;"> 55 </td>
   <td style="text-align:right;"> 10 </td>
   <td style="text-align:right;"> 81 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 29 </td>
   <td style="text-align:right;"> 82 </td>
   <td style="text-align:right;"> 21 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 59 </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 30 </td>
   <td style="text-align:right;"> 61 </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 58 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 24 </td>
   <td style="text-align:right;"> 47 </td>
   <td style="text-align:right;"> 53 </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 51 </td>
   <td style="text-align:right;"> 35 </td>
   <td style="text-align:right;"> 18 </td>
   <td style="text-align:right;"> 41 </td>
   <td style="text-align:right;"> 30 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 35 </td>
   <td style="text-align:right;"> 39 </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 49 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 79 </td>
   <td style="text-align:right;"> 33 </td>
   <td style="text-align:right;"> 62 </td>
   <td style="text-align:right;"> 45 </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 9 </td>
   <td style="text-align:right;"> 72 </td>
   <td style="text-align:right;"> 61 </td>
   <td style="text-align:right;"> 35 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 74 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 56 </td>
   <td style="text-align:right;"> 82 </td>
   <td style="text-align:right;"> 25 </td>
   <td style="text-align:right;"> 65 </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 29 </td>
   <td style="text-align:right;"> 55 </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 39 </td>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 93 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 38 </td>
   <td style="text-align:right;"> 80 </td>
   <td style="text-align:right;"> 19 </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 34 </td>
   <td style="text-align:right;"> 54 </td>
   <td style="text-align:right;"> 10 </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 72 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 40 </td>
   <td style="text-align:right;"> 80 </td>
   <td style="text-align:right;"> 106 </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 105 </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 81 </td>
   <td style="text-align:right;"> 58 </td>
   <td style="text-align:right;"> 49 </td>
   <td style="text-align:right;"> 74 </td>
   <td style="text-align:right;"> 93 </td>
   <td style="text-align:right;"> 72 </td>
   <td style="text-align:right;"> 0 </td>
  </tr>
</tbody>
</table></div>

Расстояния по автодорогам не совсем "евклидовы": например, из Кронштадта до Санкт-Петербурга нельзя добраться по прямой.

---

## Обычная проекция исходного пространства на плоскость неудачна

Карта<br/><br/>

![](02_nMDS_files/figure-html/gg-spb-1.png)
]

Проекция при помощи анализа главных координат (PCoA, многомерное шкалирование)

![](02_nMDS_files/figure-html/pcoa-spb-1.png)
]

- Карта неудачна: Стрельна, Красное Село и Гатчина сливаются вместе, хотя они не рядом.

---

## Можно восстановить карту, сохраняя ранги расстояний между объектами

Карта

![](02_nMDS_files/figure-html/gg-spb-1.png)
]

Ординация nMDS

![](02_nMDS_files/figure-html/ord-spb-1.png)
]

- Что странного в этой карте?

---

## Важные свойства nMDS-ординации

1. __Сохраняется ранг расстояний между объектами__ (похожие располагаются близко, непохожие --- далеко.  
Если объект А похож на B, больше чем на C, то и на ординации он, скорее всего, окажется ближе к B, чем к C).

2. __Имеет смысл только взаиморасположение объектов друг относительно друга__.  
Суть ординации не изменится, если решение 
  - вращать, 
  - перемещать, 
  - зеркально отражать

3. __Численные значения координат точек в ординации лишены смысла__.  
Их вообще можно не приводить на итоговой ординации.

---

# Как работает nMDS

---

## Наблюдаемые различия и ожидаемые расстояния

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

<div style="border: 1px solid #ddd; padding: 0px; overflow-y: scroll; height:450px; overflow-x: scroll; width:100%; "><table class="table" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;position: sticky; top:0; background-color: #FFFFFF;"> item1 </th>
   <th style="text-align:left;position: sticky; top:0; background-color: #FFFFFF;"> item2 </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> observed </th>
   <th style="text-align:right;position: sticky; top:0; background-color: #FFFFFF;"> nmds </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 7.090 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 9 </td>
   <td style="text-align:right;"> 9.175 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 10 </td>
   <td style="text-align:right;"> 13.668 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 18 </td>
   <td style="text-align:right;"> 20.974 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 19 </td>
   <td style="text-align:right;"> 11.966 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:right;"> 20 </td>
   <td style="text-align:right;"> 21.466 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 21 </td>
   <td style="text-align:right;"> 14.146 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 20.208 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 24 </td>
   <td style="text-align:right;"> 27.296 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:right;"> 24 </td>
   <td style="text-align:right;"> 18.554 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 25 </td>
   <td style="text-align:right;"> 23.498 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 25.232 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 26.409 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 28 </td>
   <td style="text-align:right;"> 25.045 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 29 </td>
   <td style="text-align:right;"> 24.245 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 29 </td>
   <td style="text-align:right;"> 26.272 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:right;"> 30 </td>
   <td style="text-align:right;"> 21.161 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 30 </td>
   <td style="text-align:right;"> 28.383 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 26.909 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 33.184 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 31 </td>
   <td style="text-align:right;"> 35.053 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 33 </td>
   <td style="text-align:right;"> 32.226 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 34 </td>
   <td style="text-align:right;"> 30.325 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 35 </td>
   <td style="text-align:right;"> 34.448 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 35 </td>
   <td style="text-align:right;"> 29.422 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 34.384 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 38.723 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 38.704 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 38 </td>
   <td style="text-align:right;"> 30.554 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 39 </td>
   <td style="text-align:right;"> 39.270 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 40 </td>
   <td style="text-align:right;"> 36.712 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 41 </td>
   <td style="text-align:right;"> 44.955 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:right;"> 42 </td>
   <td style="text-align:right;"> 46.194 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 43 </td>
   <td style="text-align:right;"> 43.776 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:right;"> 44 </td>
   <td style="text-align:right;"> 41.266 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 45 </td>
   <td style="text-align:right;"> 43.818 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 47 </td>
   <td style="text-align:right;"> 40.291 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 47 </td>
   <td style="text-align:right;"> 51.697 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 49 </td>
   <td style="text-align:right;"> 53.681 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 51 </td>
   <td style="text-align:right;"> 48.395 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 46.932 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 44.403 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 52 </td>
   <td style="text-align:right;"> 44.754 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Санкт-Петербург </td>
   <td style="text-align:right;"> 53 </td>
   <td style="text-align:right;"> 59.501 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:right;"> 53 </td>
   <td style="text-align:right;"> 52.822 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 54 </td>
   <td style="text-align:right;"> 51.840 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 55 </td>
   <td style="text-align:right;"> 55.352 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 56 </td>
   <td style="text-align:right;"> 53.181 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:right;"> 56 </td>
   <td style="text-align:right;"> 54.956 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:right;"> 57 </td>
   <td style="text-align:right;"> 50.175 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 58 </td>
   <td style="text-align:right;"> 54.281 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 59 </td>
   <td style="text-align:right;"> 49.355 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 61 </td>
   <td style="text-align:right;"> 57.610 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 62 </td>
   <td style="text-align:right;"> 59.808 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:right;"> 63 </td>
   <td style="text-align:right;"> 55.360 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 56.854 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 61.406 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 52.487 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 63.779 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 64 </td>
   <td style="text-align:right;"> 57.401 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 65 </td>
   <td style="text-align:right;"> 61.348 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 63.852 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 62.551 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 66 </td>
   <td style="text-align:right;"> 55.089 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 68 </td>
   <td style="text-align:right;"> 58.396 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 70.989 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Парголово </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 72.391 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 68.271 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 63.028 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 70 </td>
   <td style="text-align:right;"> 70.737 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:right;"> 71 </td>
   <td style="text-align:right;"> 70.925 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:right;"> 71 </td>
   <td style="text-align:right;"> 65.643 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 72 </td>
   <td style="text-align:right;"> 62.256 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 72 </td>
   <td style="text-align:right;"> 77.747 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Всеволожск </td>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 73 </td>
   <td style="text-align:right;"> 72.142 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 74 </td>
   <td style="text-align:right;"> 70.753 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:right;"> 75 </td>
   <td style="text-align:right;"> 69.970 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:right;"> 76 </td>
   <td style="text-align:right;"> 75.287 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:right;"> 77 </td>
   <td style="text-align:right;"> 78.277 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Сертолово </td>
   <td style="text-align:right;"> 79 </td>
   <td style="text-align:right;"> 80.049 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 80 </td>
   <td style="text-align:right;"> 71.814 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Стрельна </td>
   <td style="text-align:right;"> 80 </td>
   <td style="text-align:right;"> 78.585 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Красное Село </td>
   <td style="text-align:right;"> 81 </td>
   <td style="text-align:right;"> 84.628 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Петергоф </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 81 </td>
   <td style="text-align:right;"> 91.044 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:right;"> 82 </td>
   <td style="text-align:right;"> 82.785 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Пушкин </td>
   <td style="text-align:right;"> 82 </td>
   <td style="text-align:right;"> 86.022 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Павловск </td>
   <td style="text-align:right;"> 87 </td>
   <td style="text-align:right;"> 92.807 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Сестрорецк </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 93 </td>
   <td style="text-align:right;"> 88.408 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Гатчина </td>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:right;"> 97 </td>
   <td style="text-align:right;"> 106.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Кронштадт </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 105 </td>
   <td style="text-align:right;"> 102.031 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Зеленогорск </td>
   <td style="text-align:left;"> Шлиссельбург </td>
   <td style="text-align:right;"> 106 </td>
   <td style="text-align:right;"> 102.343 </td>
  </tr>
</tbody>
</table></div>

---

## Диаграмма Шепарда

__Диаграмма Шепарда__ --- график зависимости значений коэффициентов различия в исходном пространстве `$D_{h, i}$` и соответствующих им расстояний между точками на полученной ординации `$d_{h, i}$`.

![](02_nMDS_files/figure-html/fig-shep-1.png)

]

- По оси X --- значения коэффициентов различий `$D_{h,i}$`, ранжированные в порядке возрастания
- По оси Y --- значения расстояний на плоскости ординации --- `$d_{h,i}$`
]

В идеале, ранги расстояний между точками  `$d_{h,i}$` на хорошей ординации nMDS должны совпадать с рангами коэффициентов различия `$D_{h,i}$` в исходном многомерном пространстве.

---

## Интерпретация диаграммы Шепарда

Форма и расположение облака точек на диаграмме Шепарда позволяет судить о качестве ординации.

![](./images/shepard-plot-interpretation.png)
]

- A --- ординация объясняет большую часть исходной изменчивости
- В --- ординация объясняет небольшую часть изменчивости, но относительное расположение точек сохранено
- С --- относительное расположение точек плохо соответствует исходному

]

---

## Диаграмма Шепарда

В идеале, большим расстояниям в исходном пространстве должны соответствовать большие расстояния на ординации, а меньшим --- меньшие (т.е. должны сохраняться ранги расстояний).

![](02_nMDS_files/figure-html/fig-stress-1.png)

]

Красная линия --- предсказанные монотонной регрессией значения расстояний на ординации --- монотонно возрастает.

]

---

## Монотонная регрессия

__Монотонная регрессия__ минимизирует сумму квадратов отклонений значений координат на ординации от значений в исходном многомерном пространстве, сохраняя монотонные отношения между ними.

Алгоритм подбора монотонной регрессии называется   Pool-Adjacent-Violator-Algorithm (PAVA)

![](02_nMDS_files/figure-html/isoreg-ani-.gif)

---

## Стресс --- мера качества ординации

__Стресс (Stress)__ показывает, насколько соответствуют друг другу взаиморасположение точек на плоскости ординации и в исходном многомерном пространстве признаков.

Стресс 1 (Stress(1), Kruskal, 1964a) --- это корень из суммы квадратов остатков монотонной регрессии, отнесенной к сумме квадратов всех коэффициентов различия в исходной матрице.

![](02_nMDS_files/figure-html/fig-stress-plot-1.png)

]

`$$Stress(1) = \sqrt{\frac{\sum_{h < i}{(d_{h,i} - \hat d_{h,i})^2}}{\sum_{h < i}{d_{{h,i}}^2}}}$$`

- `$d$` --- наблюдаемое значение коэффициента различия между двумя точками;
- `$\hat d$` --- значение, предсказанное монотонной регрессией.

]

---

## Классификация стресса

Хорошо, когда `$Stress(1) < 0.2$`.

Kruskal в своей работе 1964 приводит следующую классификацию значений стресса:

Значение стресса | Оценка качества
--- | ---
0.2 | poor 
0.1 | fair
0.05 | good
0.025 | excellent
0 | "perfect"

Значение стресса в ординации с большим количеством измерений будет всегда меньше!
---

## Стресс --- не абсолютное значение

Что делать, если значение стресса в ординации в двух измерениях ~ 0.2? В работе Dexter et al., 2018 предлагают использовать пермутационный метод для оценки значения стресса в полученной ординации.

1. Выдвигается нулевая гипотеза об отсутствии какого-либо структурного распределения в ординации.

2. Значения в исходных данных меняют местами, вычисляют новую ординацию --- например, 1000 раз (пермутационный метод!).

3. Строится распределение полученных значений стресса, которое можно сравнить с найденным для исходных данных (с помощью t-теста).

---

## Стресс --- не абсолютное значение

Оценку теста с помощью пермутации можно сделать с помощью функции `oecosimu` из пакета `vegan`.

```r
oeco_res <- oecosimu(D, metaMDS, "quasiswap_count", 1000, 
                     statistic = "stress")
```

```r
hist(oeco_res$oecosimu$simulated)
abline(v = spb_ord$stress, col = "red")
```

![](02_nMDS_files/figure-html/unnamed-chunk-5-1.png)

---

## Алгоритм nMDS

<ol>
<li> Выбираем _a priori_ число измерений _k_ для финального решения (число осей).

<li> Вычисляем матрицу коэффициентов различия между объектами.

<li> Размещаем объекты в пространстве _k_ осей в изначальной конфигурации (вычисляем предсказанные значения `$d_{hi}$`). **NB!** Выбор начальной конфигурации может влиять на конечный результат: может быть найден локальный минимум значения стресса вместо общего минимума.
</ol>

- запуск анализа со случайных точек несколько раз

- запуск анализа с ординации, полученной другим методом (например, анализа главных координат)

<ol start=4>
<li> Вычисляем матрицу расстояний между точками в пространстве ординации, сравниваем полученные расстояния с исходными.

<li> Вычисляем стресс.

<li> Сдвигаем точки, чтобы минимизировать стресс (метод _steepest descent_).

<li> Повторяем шаги 3-6 многократно, чтобы избежать локальных минимумов.

<li> Обычно финальную ординацию поворачивают при помощи PCA, чтобы вдоль оси X было максимальное варьирование.
</ol>

---

## Преимущества и ограничения nMDS

### Преимущества

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

### Ограничения

- nMDS --- метод визуализации данных, а не метод анализа (см. напр. perMANOVA)
- Нужно осознанно выбирать способ трансформации данных и коэффициент различия исходя из того, какие именно свойства отношений между объектами должны быть визуализированы
- Выбранное число измерений решения влияет на результат
]

---

# Пример: Симбионты мидий

---

## Пример: Симбионты мидий

Данные можно скачать [с сайта Pangaea](https://doi.pangaea.de/10.1594/PANGAEA.870537?format=textfile) ([Krapivin 2017](https://doi.org/10.1594/PANGAEA.870537))
]

<a title="Andreas Trepte, CC BY-SA 2.5 &lt;https://creativecommons.org/licenses/by-sa/2.5&gt;, via Wikimedia Commons" href="https://commons.wikimedia.org/wiki/File:Blue_mussel_Mytilus_edulis.jpg"><img width="600" alt="Blue mussel Mytilus edulis" src="https://upload.wikimedia.org/wikipedia/commons/thumb/3/32/Blue_mussel_Mytilus_edulis.jpg/512px-Blue_mussel_Mytilus_edulis.jpg"></a>
.small[
Andreas Trepte, [CC BY-SA 2.5](https://creativecommons.org/licenses/by-sa/2.5), via Wikimedia Commons
]

---

## Открываем данные

```r
dat <- read.delim("data/Krapivin_2017_Medulis-symb.tab", skip = 36, 
                  stringsAsFactors = TRUE)
str(dat)
```

```
'data.frame':	548 obs. of  23 variables:
 $ Event                 : Factor w/ 3 levels "Gorelyi_Isl",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ Date.Time             : Factor w/ 1 level "2013-07": 1 1 1 1 1 1 1 1 1 1 ...
 $ Ord.No                : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Area                  : Factor w/ 3 levels "B Solovkecki Isl",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Site                  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Zone                  : Factor w/ 3 levels "down","middle",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Comment               : int  1 1 1 1 1 1 1 1 1 1 ...
 $ M..edulis.L.shell..mm.: num  31.2 41.1 46.4 40 34.5 47.6 34.3 49.6 49.8 30 ...
 $ M..edulis.age..a.     : num  1 2 2 3 4 3 2 4 3 1 ...
 $ Urastoma.spp.....     : int  26 16 19 11 14 16 10 4 11 0 ...
 $ Renicola.spp.....     : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Himasthla.spp.....    : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Metacercaria....      : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Gymnophallus.spp..... : int  0 0 0 0 0 0 0 0 4 0 ...
 $ Bucephalidae....      : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Alg....               : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Nematoda....          : int  1 1 0 0 0 0 0 0 0 0 ...
 $ Microsetella....      : int  0 0 0 0 0 2 0 0 0 0 ...
 $ Copepoda....          : int  0 0 0 0 0 1 0 0 0 0 ...
 $ Chironomidae....      : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Halacaridae....       : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Jaera.spp.....        : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Ostrac....            : int  0 0 0 0 0 0 0 0 0 0 ...
```

---

## Приводим в порядок названия переменных

```r
colnames(dat) # Было
```

```
 [1] "Event"                  "Date.Time"              "Ord.No"                
 [4] "Area"                   "Site"                   "Zone"                  
 [7] "Comment"                "M..edulis.L.shell..mm." "M..edulis.age..a."     
[10] "Urastoma.spp....."      "Renicola.spp....."      "Himasthla.spp....."    
[13] "Metacercaria...."       "Gymnophallus.spp....."  "Bucephalidae...."      
[16] "Alg...."                "Nematoda...."           "Microsetella...."      
[19] "Copepoda...."           "Chironomidae...."       "Halacaridae...."       
[22] "Jaera.spp....."         "Ostrac...."            
```

```r
colnames(dat) <- gsub("[.]", replacement = "", colnames(dat))
dat <- rename(dat, L = "MedulisLshellmm", Age = "Medulisagea")

colnames(dat) # Стало
```

```
 [1] "Event"           "DateTime"        "OrdNo"           "Area"           
 [5] "Site"            "Zone"            "Comment"         "L"              
 [9] "Age"             "Urastomaspp"     "Renicolaspp"     "Himasthlaspp"   
[13] "Metacercaria"    "Gymnophallusspp" "Bucephalidae"    "Alg"            
[17] "Nematoda"        "Microsetella"    "Copepoda"        "Chironomidae"   
[21] "Halacaridae"     "Jaeraspp"        "Ostrac"         
```

---

## Приводим в порядок данные

```r
# Делаем сайт фактором
dat$Site <- factor(dat$Site, levels = c(2, 3, 1), labels = c("G","L","S"))

# Сливаем редкие виды в одну категорию
f_remove <- c("Nematoda", "Microsetella", "Copepoda", "Chironomidae", 
              "Halacaridae", "Jaeraspp", "Ostrac")
dat$Other <- rowSums(dat[, f_remove])

# Суммарная численность симбионтов
f_sp <- c("Urastomaspp", "Renicolaspp", "Himasthlaspp", "Metacercaria", 
          "Gymnophallusspp", "Alg", "Other")
dat$Total  <- rowSums(dat[, f_sp])

# Данные для анализа
# Только мидии с симбионтами и возрастом от 3 до 8 лет 
dat <- dat[dat$Total != 0 & dat$Age %in% 3:8, ]
spec <- dat[, f_sp]                         # виды-симбионты
env <- dat[, c("Zone", "Site", "L", "Age")] # свойства мидий-хозяев
```

---

# Ординация сообществ симбионтов в мидиях

---

## Задание 1

Постройте ординацию мидий по обилию разных видов-симбионтов с использованием коэффициента различия Брея-Куртиса.
Следите, чтобы алгоритму удалось найти финальное решение. Если необходимо, настройте вызов `metaMDS()`

Вычислите стресс для получившейся ординации.

Нарисуйте простейший график ординации.

```
set.seed(31987)
ord_mus <- 
```

---

## Решение

```r
set.seed(31987)
ord_mus <- metaMDS(spec, dist = "bray")
```

```
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1319 
Run 1 stress 0.14 
Run 2 stress 0.1312 
... New best solution
... Procrustes: rmse 0.009127  max resid 0.09045 
Run 3 stress 0.1309 
... New best solution
... Procrustes: rmse 0.01298  max resid 0.08933 
Run 4 stress 0.1339 
Run 5 stress 0.1332 
Run 6 stress 0.1419 
Run 7 stress 0.1305 
... New best solution
... Procrustes: rmse 0.01042  max resid 0.09214 
Run 8 stress 0.1425 
Run 9 stress 0.1331 
Run 10 stress 0.1355 
Run 11 stress 0.1313 
Run 12 stress 0.1405 
Run 13 stress 0.1299 
... New best solution
... Procrustes: rmse 0.009042  max resid 0.07283 
Run 14 stress 0.1313 
Run 15 stress 0.1333 
Run 16 stress 0.132 
Run 17 stress 0.1418 
Run 18 stress 0.1345 
Run 19 stress 0.1322 
Run 20 stress 0.1308 
*** Best solution was not repeated -- monoMDS stopping criteria:
    10: stress ratio > sratmax
    10: scale factor of the gradient < sfgrmin
```

]

Все ли в порядке с этими результатами?

---

## Решение

```r
set.seed(31987)
ord_mus <- metaMDS(spec, dist = "bray")
```

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

]

```
* ## Square root transformation
* ## Wisconsin double standardization
  ## Run 0 stress 0.132 
  ## Run 1 stress 0.151 
  ## Run 2 stress 0.13 
  ## ... New best solution
  ## ... Procrustes: rmse 0.0161  max resid 0.0895 
  ## Run 3 stress 0.14 
  ...
```

Эта модель не сходится!

]

```
  ...
  ## Run 19 stress 0.133 
  ## Run 20 stress 0.138 
* ## *** No convergence -- monoMDS stopping criteria:
  ##     12: stress ratio > sratmax
  ##      8: scale factor of the gradient < sfgrmin
```

---

## Решение

После отключения автоматического подбора трансформации данных и увеличения количества попыток удается найти решение.

```r
set.seed(31987)
ord_mus <- metaMDS(spec, dist = "bray", autotransform = FALSE)
```

```
Run 0 stress 0.1412 
Run 1 stress 0.1537 
Run 2 stress 0.1505 
Run 3 stress 0.1465 
Run 4 stress 0.1412 
... New best solution
... Procrustes: rmse 0.001495  max resid 0.0125 
Run 5 stress 0.1413 
... Procrustes: rmse 0.002671  max resid 0.02478 
Run 6 stress 0.1533 
Run 7 stress 0.1436 
Run 8 stress 0.1435 
Run 9 stress 0.1412 
... Procrustes: rmse 0.0002458  max resid 0.00481 
... Similar to previous best
Run 10 stress 0.1455 
Run 11 stress 0.1413 
... Procrustes: rmse 0.001474  max resid 0.02106 
Run 12 stress 0.1479 
Run 13 stress 0.1437 
Run 14 stress 0.1447 
Run 15 stress 0.1463 
Run 16 stress 0.1434 
Run 17 stress 0.1445 
Run 18 stress 0.1447 
Run 19 stress 0.1412 
... Procrustes: rmse 0.001388  max resid 0.01208 
Run 20 stress 0.1413 
... Procrustes: rmse 0.003722  max resid 0.03408 
*** Best solution repeated 1 times
```

]

---

## Решение

```r
ord_mus$stress
```

```
[1] 0.1412
```

```r
ordiplot(ord_mus)
```

![](02_nMDS_files/figure-html/fig-mus-raw, fig-mus-bw-1.png)

Но график нужно украсить --- добавить информацию о свойствах мидий.

---

## Украшенный график

Кодируем свойства мидий: горизонт литорали цветом, а сайт --- формой.

```r
# Палитры
pal_col <- c("red", "green", "steelblue")
pal_sh <- c(1, 2, 0)

# Украшенный график
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
```

![](02_nMDS_files/figure-html/fig-mus-col-1.png)

А где же легенда?

---

## Украшенный график с легендой

```r
# Палитры
pal_col <- c("red", "green", "steelblue")
pal_sh <- c(1, 2, 0)
# Украшенный график
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])

# Легенда (пример относительного и абсолютного расположения)
legend("topleft", bty = "n",
       title = "Intertidal levels: ",
       legend = levels(env$Zone), col = pal_col, pch = 15)
legend(x = 3, y = 1.5, xjust = 1, yjust = 1,
       title = "Sites: ",
       legend = levels(env$Site), col = "black", pch = pal_sh)
```

![](02_nMDS_files/figure-html/fig-mus-col-leg-1.png)

А как понять, где какие виды паразитов?

---

## Украшенный график с центроидами видов

```r
op <- par(mar = c(3, 3, 0.1, 0.1), mgp = c(2, 1, 0))
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
legend("topleft", bty = "n",
       title = "Intertidal levels: ",
       legend = levels(env$Zone), col = pal_col, pch = 15)
legend("topright", bty = "n",
       title = "Sites: ",
       legend = levels(env$Site), col = "black", pch = pal_sh)
text(ord_mus, display = "spec", cex = 0.9, col = "grey20")
par(op)
```

![](02_nMDS_files/figure-html/fig-mus-col-leg-centr-1.png)

---

# Визуализация ординации в ggplot2

---

## Задание 2

Используя данные, приведенные ниже, постройте график ординации при помощи `ggplot2`.

Покажите 
  - точками --- мидий, 
  - текстом --- центроиды для видов-симбионтов,
  - цветом --- зону литорали, 
  - формой --- сайт,
  - размером маркеров --- размер мидий.

```r
# Координаты точек (мидий)
ord_mus_pt <- data.frame(env, scores(ord_mus, display = "sites"))
head(ord_mus_pt, 2)
```

```
  Zone Site    L Age  NMDS1   NMDS2
4 down    S 40.0   3 -1.090 -0.2404
5 down    S 34.5   4 -1.104 -0.3570
```

```r
# Координаты центроидов переменных (видов-симбионтов)
ord_mus_sp <- data.frame(scores(ord_mus, display = "species"))
ord_mus_sp$Species <- rownames(ord_mus_sp)
head(ord_mus_sp, 2)
```

```
              NMDS1   NMDS2     Species
Urastomaspp -0.9273 -0.3062 Urastomaspp
Renicolaspp  0.8644  0.9727 Renicolaspp
```

---

## Решение: График с точками

```r
gg_ord_mus <- ggplot() +
  geom_point(data = ord_mus_pt, 
             aes(x = NMDS1, y = NMDS2, colour = Zone, 
                 shape = Site, size = L), alpha = 0.5)
gg_ord_mus
```

![](02_nMDS_files/figure-html/fig-gg-ord-mus-1.png)

---

## Решение: График с центроидами видов

```r
gg_ord_mus_sp <- gg_ord_mus +
  geom_text(data = ord_mus_sp, 
            aes(x = NMDS1, y = NMDS2, 
                label = Species))
gg_ord_mus_sp
```

![](02_nMDS_files/figure-html/fig-gg-ord-mus-sp-1.png)

---

# Интерпретация ординации: envfit

---

## Анализ связи с внешними переменными c помощью функции `envfit()`

Функция `envfit()` строит регрессионную модель зависимости переменной среды от координат точек на плоскости ординации:

`$$y_i = \beta_1 \text{NMDS1}_i + \beta_2 \text{NMDS2}_i + \varepsilon_i$$`

Пермутационный тест значимости позволяет работать, когда не выполнены условия применимости регрессии.

---

## Можно подобрать сразу несколько моделей --- <br/>от каждого из факторов своя модель

```r
ef <- envfit(ord_mus, env[, c("Zone", "Site", "L", "Age")])
ef
```

```

***VECTORS

NMDS1  NMDS2   r2 Pr(>r)    
L   -1.000 -0.024 0.19  0.001 ***
Age  0.774 -0.633 0.04  0.003 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Permutation: free
Number of permutations: 999

***FACTORS:

Centroids:
           NMDS1 NMDS2
Zonedown   -0.65 -0.09
Zonemiddle -0.14 -0.05
Zoneup      0.83  0.14
SiteG       0.19  0.17
SiteL       0.23 -0.41
SiteS      -0.42  0.23

Goodness of fit:
       r2 Pr(>r)    
Zone 0.37  0.001 ***
Site 0.16  0.001 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Permutation: free
Number of permutations: 999
```

]

---

## Интерпретация результатов `envfit()` для непрерывных переменных

```r
ef$vectors
```

```
     NMDS1  NMDS2   r2 Pr(>r)    
L   -1.000 -0.024 0.19  0.001 ***
Age  0.774 -0.633 0.04  0.003 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Permutation: free
Number of permutations: 999
```

<br/>

- `MDS1` и `MDS2` --- направление вектора непрерывной переменной (длина вектора переменной пропорциональна корню из коэффициента детерминации для модели, т.е. `$R$`).  
- `r2` --- `$R^2$` коэффициент детерминации.     
- `Pr(>r)` --- пермутационная оценка статистической значимости.

---

## Интерпретация результатов `envfit()` для дискретных переменных

```r
ef$factors
```

```
Centroids:
           NMDS1 NMDS2
Zonedown   -0.65 -0.09
Zonemiddle -0.14 -0.05
Zoneup      0.83  0.14
SiteG       0.19  0.17
SiteL       0.23 -0.41
SiteS      -0.42  0.23

Goodness of fit:
       r2 Pr(>r)    
Zone 0.37  0.001 ***
Site 0.16  0.001 ***
---
Signif. codes:  
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Permutation: free
Number of permutations: 999
```

]

.pull-right[
- `MDS1` и `MDS2` --- координаты центроидов для дискретных факторов  (т.е. положение центра облака точек, относящихся к категории фактора).
- `r2` --- `$R^2$` коэффициент детерминации.    
- `Pr(>r)` --- пермутационная оценка статистической значимости. 
]

---

## График с векторами и центроидами, найденными `envfit()`

```r
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
plot(ef)
```

![](02_nMDS_files/figure-html/fig-ord-mus-envfit-1.png)

---

## Для ggplot-графика удобно использовать вспомогательный пакет

```r
# install.packages("devtools")
# devtools::install_github("gavinsimpson/ggvegan")
library(ggvegan)
ord_mus_ef <- fortify(ef)
ord_mus_ef
```

```
       label     type   NMDS1    NMDS2
1          L   Vector -0.4323 -0.01024
2        Age   Vector  0.1525 -0.12478
3   Zonedown Centroid -0.6485 -0.08635
4 Zonemiddle Centroid -0.1440 -0.05243
5     Zoneup Centroid  0.8299  0.14001
6      SiteG Centroid  0.1892  0.16946
7      SiteL Centroid  0.2281 -0.40615
8      SiteS Centroid -0.4196  0.22532
```

---

## ggplot2 версия графика `envfit()`

```r
gg_ord_mus +
  geom_segment(data = ord_mus_ef[ord_mus_ef$type == "Vector", ],
               aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
               arrow = arrow(length = unit(0.25, "cm"))) +
  geom_text(data = ord_mus_ef[ord_mus_ef$type == "Vector", ],
            aes(x = NMDS1, y = NMDS2, label = label, hjust = 1.1, vjust = 1)) +
  geom_text(data = ord_mus_ef[ord_mus_ef$type == "Centroid", ],
            aes(x = NMDS1, y = NMDS2, label = label, hjust = 1.1, vjust = 1))
```

![](02_nMDS_files/figure-html/gg-envfit-sarter-1.png)

---

## Ограничения `envfit()`

__Осторожно!__ `envfit()` основан на предположении о линейном характере связи между внешней переменной и координатами на ординации. Это не всегда так.

]
---

# Интерпретация ординации: ordisurf

---

## Анализ связи с внешними переменными c помощью `ordisurf()`

`ordisurf()` применяется, когда зависимость нелинейна.

В основе работы функции --- общая аддитивная модель (General Additive Model, GAM).

`$$Y_i = f(MDS1_i, MDS2_i)$$`

GAM позволяют подобрать нелинейные зависимости, даже если заранее неизвестна их форма.

`ordisurf()` использует пакет `mgcv` для подбора GAM на основе метода тонких пластин (thin plate spline)

---

## График с поверхностью, найденной `ordisurf()`

Изолинии на графике ординации показывают предсказанное значение зависимой переменной для разных участков ординации.

```r
par(mfrow = c(1, 2))
os_L <- ordisurf(ord_mus, env$L, method = "REML")
os_Age <- ordisurf(ord_mus, env$Age, method = "REML")
par(mfrow = c(1, 1))
```

![](02_nMDS_files/figure-html/fig-ordisurf-1.png)

---

## Интерпретация результатов `ordisurf()`

Длина мидий линейно меняется на ординации (эффективное число степеней свободы для сплайна `$edf=1.96$`). Зависимость умеренной силы (объясненная девианса 18.7%).

```r
summary(os_L)
```

```

Family: gaussian 
Link function: identity

Formula:
y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   43.410      0.484    89.7   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
          edf Ref.df   F p-value    
s(x1,x2) 1.96      9 9.8  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.183   Deviance explained = 18.7%
-REML = 1457.7  Scale est. = 92.494    n = 395
```

---

## Интерпретация результатов `ordisurf()`

Возраст мидий меняется нелинейно вдоль ординационных осей (эффективное число степеней свободы для сплайна `$edf=4.99$`). Эта зависимость довольно слабая (объясненная девианса 8.97%).

```r
summary(os_Age)
```

```

Family: gaussian 
Link function: identity

Formula:
y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.8861     0.0714    68.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
          edf Ref.df    F    p-value    
s(x1,x2) 5.02      9 3.73 0.00000099 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.0786   Deviance explained = 9.03%
-REML = 704.93  Scale est. = 2.011     n = 395
```

---

## Для ggplot-графика понадобится добыть данные о контурах

```r
fortify_ordisurf <- function(model) {
  # Fortifies an object of class ordisurf to produce 
  # a dataframe of contour coordinates 
  # for subsequent plotting in ggplot.
  xy <- expand.grid(x = model$grid$x, y = model$grid$y)
  xyz <- cbind(xy, c(model$grid$z))
  names(xyz) <- c("x", "y", "z")
  return(na.omit(xyz))
}

ord_mus_os <- fortify_ordisurf(os_Age)
head(ord_mus_os, 4)
```

```
         x      y     z
9  -0.4130 -1.371 4.826
10 -0.3089 -1.371 4.840
11 -0.2049 -1.371 4.858
12 -0.1008 -1.371 4.879
```

---

## ggplot2 версия графика с поверхностью, найденной `ordisurf()`

```r
ggplot(data = ord_mus_os, aes(x = x, y = y, z = z)) +
  stat_contour(aes(colour = after_stat(level))) +
   scale_colour_gradient(low = "blue", high = "red") +
  labs(x = "NMDS1", y = "NMDS2", colour = "Age")
```

![](02_nMDS_files/figure-html/gg-envfit-contour-1.png)

---

## Ограничения `ordisurf()`

.content-box-red[
__Осторожно!__ Для использования `ordisurf()` важно, чтобы в ординации участвовало много объектов и они располагались на плоскости ординации более-менее равномерно. Пустоты и отскоки снижают качество предсказаний.
]

---

## Финальный график

![](02_nMDS_files/figure-html/gg-envfit-fin-1.png)

---

## Take-home messages

+ nMDS --- один из методов ординации в пространстве сниженной размерности. При ординации этим методом сохраняются ранги расстояний между объектами.
+ Мера качества ординации называется "стресс".
+ Значения координат не имеют особенного смысла. Имеет значение только взаиморасположение точек.
+ Результат nMDS зависит от выбора меры различия.
+ Существует несколько способов визуализировать на ординации значения внешних переменных. При выборе между `envfit()` и `ordisurf()` обращайте внимание на форму связи.

---

## Литература

+ Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.
+ Legendre, P., Legendre, L., 2012. Numerical ecology. Elsevier.
+ Oksanen, J., 2015. Multivariate analysis of ecological communities in R: vegan tutorial. http://www.cc.oulu.fi/~jarioksa/opetus/metodi/vegantutor.pdf
+ Zuur, A. F., Ieno, E. N., Smith, G. M. Analysing Ecological Data. Springer 2007

---

# Самостоятельная работа

---

## Задание 3

**Данные**

1. Растительные сообщества во франции La Mafragh (данные из работы de Belair et al. 1987, данные `mafragh`, пакет `ade4`).
2. Деревья на острове Barro Colorado (данные из работы Condit et al. (2002), данные `BCI`, пакет `vegan`).

Во всех примерах необходимо:

1. Разобраться с данными.   
2. Построить ординацию объектов (описаний, проб и т.п.).      
3. Визуализировать связь между полученной ординацией и параметрами среды.    
4. Сделать выводы о наиболее важных факторах.