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

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

---

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

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

- Построить диаграмму 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)
]

---

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

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

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

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

nMDS ординация старается сохранить отношения между объектами.

---

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

![](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.09 </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.18 </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.67 </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.97 </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.97 </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.47 </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.15 </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.21 </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.30 </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.55 </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.50 </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.23 </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.41 </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.05 </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.25 </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.27 </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.16 </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.38 </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.91 </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.18 </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.05 </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.23 </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.32 </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.45 </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.42 </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.38 </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.72 </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.70 </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.55 </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.27 </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.71 </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.95 </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.19 </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.78 </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.27 </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.82 </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.29 </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.70 </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.68 </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.40 </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.93 </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.40 </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.75 </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.50 </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.82 </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.84 </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.35 </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.18 </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.96 </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.17 </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.28 </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.35 </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.61 </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.81 </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.36 </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.85 </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.41 </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.49 </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.78 </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.40 </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.35 </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.85 </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.55 </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.09 </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.40 </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.99 </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.39 </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.27 </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.03 </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.74 </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.92 </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.64 </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.26 </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.75 </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.14 </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.75 </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.97 </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.29 </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.28 </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.05 </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.81 </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.58 </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.63 </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.04 </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.78 </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.02 </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.81 </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.41 </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.00 </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.03 </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.34 </td>
  </tr>
</tbody>
</table></div>

---

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

__Диаграмма Шепарда__ --- график зависимости расстояний между точками на ординации и соответствующих им значений коэффициентов различия в исходном пространстве.

![](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/unnamed-chunk-3-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/unnamed-chunk-4-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$`

]

---

## Алгоритм nMDS

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

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

3. Вначале размещаем объекты в пространстве _k_ осей случайно (или PCoA).

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

5. Вычисляем стресс.

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

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

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

.small[
Марка СССР Ёжик в тумане (1988, ЦФА №5919). [Википедия](https://ru.wikipedia.org/wiki/%D0%A4%D0%B0%D0%B9%D0%BB:1988_CPA_5919.jpg)
]

---

## Преимущества и ограничения 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()`

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

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

```
ord_mus <- 
```

---

## Решение

```r
ord_mus <- metaMDS(spec, dist = "bray")
```

```
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.132 
## Run 1 stress 0.132 
## ... New best solution
## ... Procrustes: rmse 0.00977  max resid 0.09 
## Run 2 stress 0.139 
## Run 3 stress 0.139 
## Run 4 stress 0.131 
## ... New best solution
## ... Procrustes: rmse 0.0128  max resid 0.0892 
## Run 5 stress 0.143 
## Run 6 stress 0.142 
## Run 7 stress 0.134 
## Run 8 stress 0.131 
## ... New best solution
## ... Procrustes: rmse 0.0135  max resid 0.0895 
## Run 9 stress 0.134 
## Run 10 stress 0.131 
## ... New best solution
## ... Procrustes: rmse 0.00869  max resid 0.0824 
## Run 11 stress 0.131 
## ... New best solution
## ... Procrustes: rmse 0.00772  max resid 0.0753 
## Run 12 stress 0.137 
## Run 13 stress 0.13 
## ... New best solution
## ... Procrustes: rmse 0.0114  max resid 0.0884 
## Run 14 stress 0.141 
## Run 15 stress 0.131 
## Run 16 stress 0.132 
## Run 17 stress 0.134 
## Run 18 stress 0.142 
## Run 19 stress 0.131 
## Run 20 stress 0.132 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     10: stress ratio > sratmax
##      9: scale factor of the gradient < sfgrmin
```

]

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

---

## Решение

```r
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
ord_mus <- metaMDS(spec, dist = "bray", autotransform = FALSE)
```

```
## Run 0 stress 0.141 
## Run 1 stress 0.141 
## ... Procrustes: rmse 0.00177  max resid 0.0334 
## Run 2 stress 0.141 
## ... Procrustes: rmse 0.00254  max resid 0.0252 
## Run 3 stress 0.141 
## ... Procrustes: rmse 0.0029  max resid 0.0337 
## Run 4 stress 0.146 
## Run 5 stress 0.141 
## ... New best solution
## ... Procrustes: rmse 0.00136  max resid 0.0132 
## Run 6 stress 0.144 
## Run 7 stress 0.142 
## ... Procrustes: rmse 0.00489  max resid 0.0719 
## Run 8 stress 0.154 
## Run 9 stress 0.146 
## Run 10 stress 0.157 
## Run 11 stress 0.141 
## ... Procrustes: rmse 0.0016  max resid 0.0157 
## Run 12 stress 0.141 
## ... Procrustes: rmse 0.000587  max resid 0.00906 
## ... Similar to previous best
## Run 13 stress 0.143 
## Run 14 stress 0.141 
## ... Procrustes: rmse 0.00178  max resid 0.0334 
## Run 15 stress 0.141 
## ... Procrustes: rmse 0.00334  max resid 0.0334 
## Run 16 stress 0.141 
## ... Procrustes: rmse 0.00227  max resid 0.0247 
## Run 17 stress 0.141 
## ... Procrustes: rmse 0.00133  max resid 0.0128 
## Run 18 stress 0.141 
## ... Procrustes: rmse 0.00128  max resid 0.0128 
## Run 19 stress 0.141 
## ... Procrustes: rmse 0.00361  max resid 0.0331 
## Run 20 stress 0.141 
## ... Procrustes: rmse 0.00187  max resid 0.0333 
## *** Solution reached
```

]

---

## Решение

```r
ord_mus$stress
```

```
## [1] 0.141
```

```r
ordiplot(ord_mus)
```

![](02_nMDS_files/figure-html/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.09 -0.240
## 5 down    S 34.5   4 -1.10 -0.357
```

```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.926 -0.306 Urastomaspp
## Renicolaspp  0.863  0.971 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.025 0.19  0.001 ***
## Age  0.776 -0.631 0.04  0.001 ***
## ---
## 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.025 0.19  0.001 ***
## Age  0.776 -0.631 0.04  0.001 ***
## ---
## 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/unnamed-chunk-15-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.433 -0.0107
## 2        Age   Vector  0.152 -0.1238
## 3   Zonedown Centroid -0.648 -0.0866
## 4 Zonemiddle Centroid -0.144 -0.0531
## 5     Zoneup Centroid  0.829  0.1409
## 6      SiteG Centroid  0.189  0.1702
## 7      SiteL Centroid  0.228 -0.4074
## 8      SiteS Centroid -0.420  0.2258
```

---

## 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/unnamed-chunk-17-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/unnamed-chunk-18-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.81  <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.6  Scale est. = 92.464    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.0713    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.03      9 3.77 0.00000089 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0793   Deviance explained = 9.11%
## -REML =  704.8  Scale est. = 2.0093    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.429 -1.37 4.84
## 10 -0.328 -1.37 4.85
## 11 -0.226 -1.37 4.87
## 12 -0.124 -1.37 4.88
```

---

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

```r
ggplot(data = ord_mus_os, aes(x = x, y = y, z = z)) +
  stat_contour(aes(colour = ..level..)) +
  labs(x = "NMDS1", y = "NMDS2", colour = "Age")
```

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

---

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

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

---

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

![](02_nMDS_files/figure-html/unnamed-chunk-23-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. Сделать выводы о наиболее важных факторах.