Анализ избыточности (Redundancy analysis, RDA)

  • Связь нескольких наборов переменных
  • Анализ избыточности, теория и практика
  • Проверка значимости ординации
  • Выбор оптимальной модели
  • Частный анализ избыточности и компоненты объясненной инерции
  • Компоненты объясненной изменчивости

Вы сможете

  • Проводить анализ избыточности
  • Оценивать долю объясненной инерции
  • Интерпретировать компоненты по нагрузкам переменных
  • Строить ординацию объектов в пространстве компонент
  • Проверять значимость модели ординации при помощи пермутационного теста
  • Разделять объясненную инерцию на компоненты, связанные с разными наборами переменных, при помощи частного анализа избыточности

Связь нескольких наборов переменных

Что будет определять генетическую структуру в колониях бабочек?

Пример: генетика бабочек Euphydryas editha

Частоты разных аллелей фосфоглюкоизомеразы и данные о факторах среды для 16 колоний бабочек Euphydryas editha в Калифорнии и Орегоне (данные McKechnie et al., 1975)

Winged Wonder by Roger Lynn on Flickr

Анализ избыточности

Анализ избыточности (RDA, Redundancy analysis)

RDA — метод прямой ординации (ограниченной ординации = constrained ordination). Предложен Rao в 1964 г., переоткрыт Wollenberg в 1977 г.

  • Прямое расширение возможностей регрессионного анализа на многомерные данные.

  • Основан на анализе главных компонент.

  • Нужно две матрицы данных: матрица зависимых переменных и матрица предикторов.

  • Нужно найти такие компоненты матрицы зависимых переменных, которые являются линейными комбинациями предикторов и отражают максимум изменчивости.

Непрямая ординация и множественная регрессия

Рис. из Legendre, Legendre, 1998

Прямая (ограниченная) ординация

Рис. из Legendre, Legendre, 1998

Если множественную линейную регрессию можно обобщить для нескольких зависимых переменных — получится RDA

\[y _{i} = b _0 + b _1 x _{1i} + b _2 x _{2i} + ... + b _k x_{ki} + \epsilon _{i}\]

Уравнение множественной линейной регрессии можно переписать в виде матриц.

\[\left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right] = \left[\begin{array}{cc} 1 & x_{1,1} & x_{1,2} & \cdots & x_{1,k} \\ 1 & x_{2,1} & x_{2,2} & \cdots & x_{2,k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n,1} & x_{n,2} & \cdots & x_{n,k} \end{array}\right] \cdot \left[\begin{array}{c} b _0 \\ b _1 \\ b _2 \\ \vdots \\ b _k \end{array}\right] + \left[\begin{array}{c} \epsilon _1 \\ \epsilon _2 \\ \vdots \\ \epsilon _n \end{array}\right] \]

Сокращенная форма записи: \(\mathbf{y} = \mathbf{X} \mathbf{b} + \mathbf{\epsilon}\), причем \(\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}(\mathbf{X}' \mathbf{y})\).

Во множественной линейной регрессии \(\mathbf{\hat y} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}(\mathbf{X}'\mathbf{y})\)

В RDA зависимая переменная — матрица, т.е. \(\mathbf{\hat Y} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}(\mathbf{X}'\mathbf{Y})\),

Подготовка данных к RDA

  • Значения откликов в матрице Y центрируют и стандартизуют при необходимости (как перед PCA).
  • Предикторы в матрице X центрируют для простоты интерпретации (т.к. в регрессии пропадут интерсепты. Иногда имеет смысл стандартизовать предикторы и/или трансформировать их.

Фрагмент Рис.11.2 из Legendre, Legendre, 1998

Вычисление канонических (constrained) осей

Фрагмент Рис.11.2 из Legendre, Legendre, 1998


1.Строим регрессию \(\mathbf{Y}_{n \times p}\) от \(\mathbf{X}_{n \times m}\), получаем предсказанные значения
\(\mathbf{\hat Y} = \mathbf{XB} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}(\mathbf{X}'\mathbf{Y})\).

2.Делаем PCA по \(\mathbf{\hat Y}\), получаем собственные векторы \(\mathbf{U}_{p \times p}\) (канонические). Всего не больше чем \(min[p, m, n-1]\) шт.

  • Ординация объектов в пространстве наблюдаемых значений: \(\mathbf{F} = \mathbf{YU}\)
  • Ординация объектов в пространстве предикторов: \(\mathbf{Z} = \mathbf{\hat YU} = XBU\)

Вычисление неканонических (unconstrained) осей

3.Делаем PCA по остаткам от множественной регрессии \(\mathbf {Y_{res}} = \mathbf{Y - \hat Y}\), получаем собственные векторы \(\mathbf{U_{res}}\) (неканонические). Всего не больше \(min[p, n - 1]\) шт.

  • Ординация объектов в пространстве остатков от регрессии: \(\mathbf{W} = \mathbf{Y_{res}U_{res}}\)

Фрагмент Рис.11.2 из Legendre, Legendre, 1998

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

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

  • Независимость наблюдений (сайтов, объектов).
  • Линейная зависимость переменных-откликов от предикторов.
  • Не должно быть мультиколлинеарности предикторов.
  • Наблюдений должно быть значительно больше, чем предикторов (как в регрессии, см. “проклятие размерности”).

Готовим к RDA данные из примера

Подготовка исходных данных

Нужно две матрицы: зависимые переменные (генетические данные для разных сайтов) и предикторы (данные об условиях среды, координатах и т.п.).

Структура исходных данных

  • $xy — координаты колоний
  • $envir — 4 фактора среды для колоний
  • $genet — частоты 6 аллелей в колониях
  • $contour — карта Калифорнии
library(ade4)
data(butterfly)
str(butterfly)
## List of 5
##  $ xy     :'data.frame': 16 obs. of  2 variables:
##   ..$ x: num [1:16] 41 57 56 57 58 67 73 165 167 77 ...
##   ..$ y: num [1:16] 238 134 131 127 124 126 121 23 18 127 ...
##  $ envir  :'data.frame': 16 obs. of  4 variables:
##   ..$ Altitude     : num [1:16] 500 800 570 550 550 380 930 650 600 1500 ...
##   ..$ Precipitation: num [1:16] 43 20 28 28 28 15 21 10 10 19 ...
##   ..$ Temp_Max     : num [1:16] 98 92 98 98 98 99 99 101 101 99 ...
##   ..$ Temp_Min     : num [1:16] 17 32 26 26 26 28 28 27 27 23 ...
##  $ genet  :'data.frame': 16 obs. of  6 variables:
##   ..$ 0.4 : num [1:16] 0 0 0 0 0 0 0 10 14 0 ...
##   ..$ 0.6 : num [1:16] 3 16 6 4 1 2 0 21 26 1 ...
##   ..$ 0.8 : num [1:16] 22 20 28 19 8 19 15 40 32 6 ...
##   ..$ 1   : num [1:16] 57 38 46 47 50 44 50 25 28 80 ...
##   ..$ 1.16: num [1:16] 17 13 17 27 35 32 27 4 0 12 ...
##   ..$ 1.3 : num [1:16] 1 13 3 3 6 3 8 0 0 1 ...
##  $ contour:'data.frame': 59 obs. of  4 variables:
##   ..$ x1: num [1:59] 29 103 102 159 211 ...
##   ..$ y1: num [1:59] 234.1 232.6 164.1 113.3 61.9 ...
##   ..$ x2: num [1:59] 103 102 159 211 210 ...
##   ..$ y2: num [1:59] 232.6 164.1 113.3 61.9 54.8 ...
##  $ Spatial:Formal class 'SpatialPolygons' [package "sp"] with 4 slots
##   .. ..@ polygons   :List of 1
##   .. .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. .. .. .. ..@ Polygons :List of 1
##   .. .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. .. .. .. ..@ labpt  : num [1:2] 110 122
##   .. .. .. .. .. .. .. ..@ area   : num 17741
##   .. .. .. .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. .. .. .. ..@ coords : num [1:60, 1:2] 29 103 102 159 211 ...
##   .. .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x1" "y1"
##   .. .. .. .. ..@ plotOrder: int 1
##   .. .. .. .. ..@ labpt    : num [1:2] 110 122
##   .. .. .. .. ..@ ID       : chr "contour"
##   .. .. .. .. ..@ area     : num 17741
##   .. ..@ plotOrder  : int 1
##   .. ..@ bbox       : num [1:2, 1:2] 23.3 13.5 211.2 234.1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "x" "y"
##   .. .. .. ..$ : chr [1:2] "min" "max"
##   .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. ..@ projargs: chr NA

Создадим переменные с более короткими названиями для удобства

gen <- butterfly$genet
head(gen, 1)
##    0.4 0.6 0.8  1 1.16 1.3
## SS   0   3  22 57   17   1
env_geo <- cbind(butterfly$envir, butterfly$xy)
head(env_geo, 1)
##    Altitude Precipitation Temp_Max Temp_Min  x   y
## SS      500            43       98       17 41 238

Нужно ли стандартизовать зависимые переменные?

summary(gen)
##       0.4             0.6             0.8             1       
##  Min.   : 0.00   Min.   : 0.00   Min.   : 1.0   Min.   :25.0  
##  1st Qu.: 0.00   1st Qu.: 2.75   1st Qu.:12.5   1st Qu.:36.8  
##  Median : 0.00   Median : 4.50   Median :18.0   Median :47.0  
##  Mean   : 1.75   Mean   : 7.19   Mean   :18.6   Mean   :51.2  
##  3rd Qu.: 0.25   3rd Qu.: 7.50   3rd Qu.:23.5   3rd Qu.:59.2  
##  Max.   :14.00   Max.   :26.00   Max.   :40.0   Max.   :92.0  
##       1.16           1.3       
##  Min.   : 0.0   Min.   : 0.00  
##  1st Qu.:10.0   1st Qu.: 0.00  
##  Median :17.0   Median : 3.00  
##  Mean   :17.2   Mean   : 4.12  
##  3rd Qu.:27.0   3rd Qu.: 6.50  
##  Max.   :35.0   Max.   :14.00
  • Нет, можно оставить так, т.к. масштаб не сильно различается

Линейны ли связи между переменными?

pairs(env_geo)

Удаляем коллинеарные предикторы

Не будем учитывать предикторы у которых VIF >= 10 (Borcard et al. 2011)

library(car)
vif(lm(gen$`0.4` ~ ., data = env_geo))
##      Altitude Precipitation      Temp_Max      Temp_Min             x 
##         23.89          4.46          4.77         22.49          7.54 
##             y 
##         11.25
vif(lm(gen$`0.4` ~ . -Temp_Min, data = env_geo))
##      Altitude Precipitation      Temp_Max             x             y 
##          7.52          4.41          4.61          5.61          7.27
  • Минимальная температура связана с высотой. Придется оставить что-то одно. Не будем использовать минимальную температуру.

RDA в R

RDA в vegan

  • Зависимые переменные (отклики) — генетические данные
  • Независимые переменные (предикторы) — переменные среды
library(vegan)
bf_rda <- rda(gen ~ Altitude + Precipitation + Temp_Max, data = env_geo)
summary(bf_rda) 
## 
## Call:
## rda(formula = gen ~ Altitude + Precipitation + Temp_Max, data = env_geo) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total             730      1.000
## Constrained       372      0.509
## Unconstrained     358      0.491
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          RDA1    RDA2    RDA3     PC1     PC2     PC3
## Eigenvalue            352.481 18.4680 0.78079 212.437 108.097 27.2173
## Proportion Explained    0.483  0.0253 0.00107   0.291   0.148  0.0373
## Cumulative Proportion   0.483  0.5084 0.50950   0.801   0.949  0.9861
##                          PC4    PC5
## Eigenvalue            8.9505 1.1684
## Proportion Explained  0.0123 0.0016
## Cumulative Proportion 0.9984 1.0000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                          RDA1    RDA2   RDA3
## Eigenvalue            352.481 18.4680 0.7808
## Proportion Explained    0.948  0.0497 0.0021
## Cumulative Proportion   0.948  0.9979 1.0000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  10.2 
## 
## 
## Species scores
## 
##        RDA1    RDA2    RDA3    PC1     PC2    PC3
## 0.4   0.679 -0.5790  0.2162 -0.947 -0.5512  0.187
## 0.6   0.910 -0.7304 -0.1210 -2.230 -0.7603  1.000
## 0.8   2.617 -0.0559 -0.1852 -2.772  0.0314 -1.414
## 1    -6.306  0.1200 -0.0531  3.396 -2.5090 -0.302
## 1.16  1.515  1.3249  0.0315  2.304  2.6504 -0.300
## 1.3   0.586 -0.0796  0.1115  0.249  1.1387  0.829
## 
## 
## Site scores (weighted sums of species scores)
## 
##        RDA1    RDA2   RDA3      PC1     PC2    PC3
## SS  -0.6880   1.729 -10.47 -0.94834 -1.0178 -2.067
## SB   1.8324  -5.038  -1.36 -2.06492  1.8111  5.099
## WSB  1.0724   0.216 -16.42 -1.08568 -0.0967 -2.968
## JRC  0.7453   6.006   2.97  0.74940  1.5211 -1.058
## JRH  0.0219  11.114  28.23  2.93439  2.7488  1.948
## SJ   1.2334   8.924   7.96  2.68930  1.8480 -2.301
## CR   0.1498   7.180  17.47  2.74983  1.1213 -0.327
## UO   4.2939 -13.727 -29.85 -4.31586 -1.9736 -1.956
## LO   3.5291 -17.677 -16.85 -4.04222 -3.5723  2.512
## DP  -4.5721   1.158   5.85  5.84888 -6.7130  0.182
## PZ   3.1672   2.245 -11.69 -1.50605  2.6764 -5.186
## MC  -2.2846  -0.769  -7.95 -1.50866 -2.2541  1.324
## IF   0.5063   1.559   7.48 -0.00751  1.0052  2.605
## AF   2.7351   3.555  26.91  0.74623  3.1376  3.036
## GH  -5.1987  -3.758  -5.28 -0.10697 -0.7903 -0.573
## GL  -6.5436  -2.717   3.01 -0.13183  0.5484 -0.270
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##       RDA1   RDA2   RDA3      PC1     PC2    PC3
## SS  -0.839  4.576 -1.590 -0.94834 -1.0178 -2.067
## SB   0.154 -0.717 -5.823 -2.06492  1.8111  5.099
## WSB  0.511  1.369 -1.733 -1.08568 -0.0967 -2.968
## JRC  0.516  1.387 -1.762  0.74940  1.5211 -1.058
## JRH  0.516  1.387 -1.762  2.93439  2.7488  1.948
## SJ   1.915 -1.155 -1.512  2.68930  1.8480 -2.301
## CR   1.232 -0.384 -0.598  2.74983  1.1213 -0.327
## UO   2.647 -2.374  0.239 -4.31586 -1.9736 -1.956
## LO   2.660 -2.330  0.165 -4.04222 -3.5723  2.512
## DP   1.275 -1.308  0.213  5.84888 -6.7130  0.182
## PZ   1.283 -0.832  2.066 -1.50605  2.6764 -5.186
## MC  -2.232  6.462  2.312 -1.50866 -2.2541  1.324
## IF   0.176  1.054  4.091 -0.00751  1.0052  2.605
## AF   1.999 -1.127  5.286  0.74623  3.1376  3.036
## GH  -4.958 -2.619 -0.750 -0.10697 -0.7903 -0.573
## GL  -6.853 -3.390  1.158 -0.13183  0.5484 -0.270
## 
## 
## Biplot scores for constraining variables
## 
##                 RDA1   RDA2  RDA3 PC1 PC2 PC3
## Altitude      -0.885 -0.403 0.231   0   0   0
## Precipitation -0.834  0.524 0.175   0   0   0
## Temp_Max       0.868  0.350 0.353   0   0   0

Структура общей изменчивости

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

## Partitioning of variance:
##               Inertia Proportion
## Total             730      1.000
## Constrained       372      0.509
## Unconstrained     358      0.491
  • Total — для всех осей, общая изменчивость исходной матрицы откликов (генетич. структуры в разных сайтах)
  • Constrained — для осей, кот. являются комбинациями предикторов (изменчивость объясненная средой)
  • Unconstrained — необъясненная предикторами изменчивость

Важность различных компонент

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

## Eigenvalues, and their contribution to the variance
## $importance
## Importance of components:
##                          RDA1    RDA2    RDA3     PC1     PC2     PC3
## Eigenvalue            352.481 18.4680 0.78079 212.437 108.097 27.2173
## Proportion Explained    0.483  0.0253 0.00107   0.291   0.148  0.0373
## Cumulative Proportion   0.483  0.5084 0.50950   0.801   0.949  0.9861
##                          PC4    PC5
## Eigenvalue            8.9505 1.1684
## Proportion Explained  0.0123 0.0016
## Cumulative Proportion 0.9984 1.0000
  • Много изменчивости объяснено, но много осталось необъясненной. Первые две канонических оси объясняют 51% изменчивости, но первые две неканонических объясняют еще 43%

Распределение изменчивости, потенциально объяснимой факторами

## Accumulated constrained eigenvalues
## $importance
## Importance of components:
##                          RDA1    RDA2   RDA3
## Eigenvalue            352.481 18.4680 0.7808
## Proportion Explained    0.948  0.0497 0.0021
## Cumulative Proportion   0.948  0.9979 1.0000
  • Первая ограниченная ось объясняет большую часть потенциально объяснимой изменчивости. Остальные оси почти ничего не объясняют.

Собственные векторы, нагрузки переменных = “species scores”

scores(bf_rda, display = "species", choices = 1:5)
##        RDA1    RDA2    RDA3    PC1     PC2
## 0.4   0.679 -0.5790  0.2162 -0.947 -0.5512
## 0.6   0.910 -0.7304 -0.1210 -2.230 -0.7603
## 0.8   2.617 -0.0559 -0.1852 -2.772  0.0314
## 1    -6.306  0.1200 -0.0531  3.396 -2.5090
## 1.16  1.515  1.3249  0.0315  2.304  2.6504
## 1.3   0.586 -0.0796  0.1115  0.249  1.1387
## attr(,"const")
## [1] 10.2

Корреляции между откликами и предикторами

  • Сильная корреляция между генетической структурой и средой только для первой ограниченной оси. Для других — умеренные или слабые.
spenvcor(bf_rda)
##  RDA1  RDA2  RDA3 
## 0.830 0.348 0.166
  • Будьте осторожны с интерпретацией! Возможна ситуация, когда корреляция между откликами и предикторами будет сильной, но при этом соотв. канонические оси будут объяснять ничтожную долю изменчивости.

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

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

  • Какие предикторы важнее всего?
  • Какими факторами определяется значение зависимых переменных?

Триплоты:

  • переменные-отклики (“species”),
  • объекты (“sites”)
  • переменные-предикторы (непрерывные в виде векторов, дискретные в виде центроидов)

Биплоты:

  • отклики + предикторы
  • объекты + предикторы

Триплот корреляций (scaling = 2): Какие переменные среды сильнее всего определяют сходство объектов?

plot(bf_rda, scaling = 2)

  • Векторы — независимые переменные, факторы среды
  • Надписи — объекты (сайты, особи, популяции и пр.)

  • Красные надписи —- зависимые переменные

  • Косинусы углов между векторами — корреляции между соотв. переменными
  • Расстояния между объектами не имеют конкретного смысла

  • Проекция объекта на линию-вектор — значение переменной для данного объекта

Пример интерпретации триплота корреляций

plot(bf_rda, scaling = 2)

  • Вдоль первой оси изменяется температура, высота и осадки
  • Вдоль второй оси — немного меняется уровень осадков

Триплот расстояний (scaling = 1): Насколько похожи друг на друга объекты?

plot(bf_rda, scaling = 1)


  • Надписи — объекты (сайты, особи, популяции и пр.)
  • Красные надписи — зависимые переменные

  • Векторы — независимые переменные, факторы среды

  • Расстояния между точками — расстояния между наблюдениями
  • Косинусы углов между векторами предикторов и откликов — корреляции между соотв. переменными, другие углы не имеют смысла
  • Проекция объекта на линию-вектор отражает примерное относительное положение данного объекта вдоль соотв. переменной (но не значение)

  • Отношения между дискретными и непрерывными предикторами не интерпретируются

Пример интерпретации триплота расстояний

plot(bf_rda, scaling = 1)

  • Генетическая структура в LO и UO похожа, но не похожа на остальные места
  • GL и GH — более высокогорные сайты, чем LO и UO

Проверка значимости ординации

Проверить значимость ординации можно несколькими способами

Каждый из способов дает ответ на свой собственный вопрос.

Влияют ли факторы на зависимые переменные?

общий тест значимости ординации

Какие именно факторы влияют на зависимые переменные?

тест факторов по отдельности (Type I или Type III эффекты)

Значима ли изменчивость вдоль осей ординации?

тест каждой из осей ординации

Псевдо-F статистика

Псевдо-F статистика оценивает соотношение изменчивости, объясненной предикторами, и остаточной изменчивости. Она рассчитывается на основе сумм собственных чисел (инерций) канонических и неканонических осей.

\[pseudoF = \frac {\Lambda_{c} / p} {\Lambda_{r}/(n - p - 1)}\]

  • \(\Lambda_{c}\) — инерция канонических осей
  • \(\Lambda_{r}\) — остаточная инерция
    тогда \(\Lambda_{c} + \Lambda_{r} = \Lambda\) — общая инерция
  • \(p\) — ранг канонических осей, \(n\) — число наблюдений.

Псевдо-F построена по тем же принципам, что F в ANOVA, но не следует F-распределению (кроме единственного частного случая RDA c одной переменной). Поэтому нам нужны пермутации, чтобы оценить форму ее распределения.

Пермутационный тест

  • Делаем случайную перестановку данных в матрице зависимых переменных
  • Рассчитываем \(pseudoF\) для пермутации.
  • Повторяем процесс заданное количество раз.

По распределению пермутационной статистики оцениваем долю пермутаций, в которых качество модели оказалось лучше, чем для реальной модели — это доверительная вероятность \(p\). Если в большинстве пермутаций получаются модели хуже, чем реальная модель, то значит зависимость от предикторов значима.

Общий тест: Влияют ли факторы на зависимые переменные?

  • тестируем гипотезу о том, что отношения между генотипом и средой значимы.

\(H _0\): значения откликов в пробах не зависят от переменных среды (генетическая структура не зависит от среды)

anova(bf_rda, permutations = 9999)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = gen ~ Altitude + Precipitation + Temp_Max, data = env_geo)
##          Df Variance    F Pr(>F)   
## Model     3      372 4.15  0.008 **
## Residual 12      358               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Есть ли связь генетической структуры со средой?

  • связь генетической структуры и среды значима

Тест факторов, type I эффекты: Какие факторы влияют на зависимые переменные?

anova(bf_rda, by = "term", permutations = 9999)
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = gen ~ Altitude + Precipitation + Temp_Max, data = env_geo)
##               Df Variance    F Pr(>F)   
## Altitude       1      279 9.37 0.0017 **
## Precipitation  1       72 2.43 0.0967 . 
## Temp_Max       1       20 0.67 0.5142   
## Residual      12      358               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Генетическая структура популяций бабочек достоверно зависит от высоты, если в модель включены др. факторы.
  • Но это Type I эффекты — они зависят от порядка включения факторов в модель. Т.е. после включения высоты в модель другие факторы уже не влияют.

Тест факторов, type III эффекты: Какие факторы влияют на зависимые переменные?

anova(bf_rda, by = "mar", permutations = 9999)
## Permutation test for rda under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = gen ~ Altitude + Precipitation + Temp_Max, data = env_geo)
##               Df Variance    F Pr(>F)
## Altitude       1       12 0.41   0.67
## Precipitation  1       71 2.39   0.10
## Temp_Max       1       20 0.67   0.51
## Residual      12      358
  • Если протестировать каждый из факторов отдельно, при условии, что все остальные включены в модель, то получится, что ни один из них не влияет.

Тест значимости осей, ограниченных факторами: Вдоль какой из осей значимо меняется генетическая структура?

  • \(H _0\): значения переменных-откликов для объектов не зависят от переменных-предикторов
  • пермутационный: выбирает оси, которые объясняют больше изменчивости, чем из др. матриц, полученных путем перестановок
anova(bf_rda, by = "axis", permutations = 9999)
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = gen ~ Altitude + Precipitation + Temp_Max, data = env_geo)
##          Df Variance     F Pr(>F)   
## RDA1      1      352 11.82 0.0054 **
## RDA2      1       18  0.62 0.8599   
## RDA3      1        1  0.03 0.9967   
## Residual 12      358                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Генетическая структура значимо меняется вдоль первой главной оси

Выбор оптимальной модели

Выбор оптимальной модели

У нас проблема. Если мы тестируем любой из факторов, после включения остальных в модель — он не влияет. Вероятно, модель не оптимальна.

Как подобрать оптимальную модель?

  • Можно использовать пошаговый выбор модели: добавляем в модель лучшие переменные и снова исключаем те, что потеряли значимость. (Вспомните, как это было для регрессионных моделей.)

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

  • Модели с разным числом предикторов можно сравнить при помощи пермутационного теста (AIC для канонических ординаций не существует!)
  • Осторожно!
  • Пошаговый выбор модели — не панацея, т.к. разные пошаговые методы могут дать разные конечные модели!
  • Не запихивайте в модель сразу все предикторы, а начинайте ваш выбор с небольшого числа важных факторов.

Пошаговый выбор оптимальной модели

Для пошагового выбора нам понадобятся полная и нулевая модели

m1 <-rda(gen ~ Altitude + Precipitation + Temp_Max, data = env_geo)
m0 <- rda(gen ~ 1, data = env_geo)

Запускаем пошаговый выбор

m <- ordistep(m0, scope = formula(m1), permutations = 9999)
## 
## Start: gen ~ 1 
## 
##                 Df AIC    F Pr(>F)   
## + Altitude       1 101 8.69 0.0023 **
## + Temp_Max       1 101 8.12 0.0031 **
## + Precipitation  1 102 7.30 0.0043 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: gen ~ Altitude 
## 
##            Df AIC    F Pr(>F)   
## - Altitude  1 106 8.69 0.0016 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                 Df   AIC    F Pr(>F)  
## + Precipitation  1  99.9 2.49  0.088 .
## + Temp_Max       1 102.0 0.64  0.551  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Оптимальная модель, отобранная при помощи пошагового алгоритма

m$anova
##            Df AIC    F Pr(>F)   
## + Altitude  1 101 8.69 0.0023 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Оптимальная модель содержит только один предиктор — высоту

Частный анализ избыточности

Зачем нужен частный анализ избыточности?

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

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

И теперь у нас два вопроса:

Кто виноват?

  • сходный климат в близких локациях
  • поток генов между близкими колониями облегчен

Что делать?

  • Нужно удалить влияние географического положения, чтобы сделать корректный вывод о связи генотипов со средой.

Частный анализ избыточности

  • зависимость от одного набора переменных (предикторов), когда влияние другого (ковариат) исключено.

Техника:

  1. Множественная регрессия зависимости предикторов от ковариат.

  2. Остатки от этой регрессии — это то, что от ковариат не зависит — можно использовать в PCA в качестве предикторов (вместо исходных переменных).

Legendre & Legendre 1998

Делаем частный RDA: зависимость генетической структуры от среды с учетом географического положения

bf_prda_1 <- rda(gen ~ Altitude + Condition(x + y), data = env_geo)
anova(bf_prda_1, permutations = 9999) ## Пермутационный тест
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = gen ~ Altitude + Condition(x + y), data = env_geo)
##          Df Variance    F Pr(>F)   
## Model     1      213 9.94 0.0014 **
## Residual 12      257               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Высота объясняет генетическую изменчивость, даже после удаления влияния географической близости

График ординации

op <- par(mfrow = c(1, 2))
plot(bf_prda_1, main = "Correlation triplot", scaling = 2)
plot(bf_prda_1, main = "Distance triplot", scaling = 1)

par(op)

Компоненты объясненной инерции

Компоненты инерции

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

Но какая часть изменчивости генетической структуры объясняется в чистом виде географической близостью, а какая — общим действием среды и географии?

Общую инерцию делим на части

Компонент инерции Что описывает? Инерция
a + b + c вся потенциально объяснимая средой и географией изменчивость Объясненная каноническими осями
RDA среда + география
a изменчивость, объясненная средой Объясненная каноническими осями
pRDA среда + Condition(география)
с изменчивость, объясненная географией Объясненная каноническими осями
pRDA география + Condition(среда)
b изменчивость, совместно объясненная средой и географией Находим инерцию по разности

Подбираем модели RDA, нужные для поиска компонентов инерции

Нам нужна полная модель RDA: генетика от среды и географического положения (чтобы найти a + b + c)

У нас уже есть частный RDA №1: зависимость генетики от среды с учетом географии (чтобы найти a)

Нам нужен частный RDA №2: генетика от географии с учетом свойств среды (чтобы найти c)

bf_prda_2 <- rda(gen ~ x + y + Condition(Altitude), data = env_geo)
bf_rda_full <- rda(gen ~ x + y + Altitude, data = env_geo)

Задание: Найдите компоненты инерции

По результатам трех RDA найдите

всю потенциально объяснимую инерцию,

а так же долю инерции, объясненную:

  • средой и географией
  • средой, но не с географией
  • географией, но не со средой

Решение: 1) Сколько инерции потенциально объясняется средой и географией?

sum_full <- summary(bf_rda_full)
> sum_full
## Partitioning of variance:
##               Inertia Proportion
## Total             730      1.000
## Constrained       473      0.648
## Unconstrained     257      0.352

Инерция, объясненная вместе средой и географией, здесь достаточно велика — 472.64

(I_total <- sum_full$constr.chi)
## [1] 473

В отличие от нее, доля инерции, объясненной ограниченной матрицей, может быть довольно малой по отношению к общей инерции, поэтому можно сосредоточиться на доле от потенциально объяснимой инерции (от sum_full$constr.chi)

Решение: 2) Инерция, объясненная средой

sum_prda_1 <- summary(bf_prda_1)
> sum_prda_1
partit(sum_prda_1)
## Partitioning of variance:
##               Inertia Proportion
## Total             730      1.000
## Conditioned       260      0.356
## Constrained       213      0.292
## Unconstrained     257      0.352
  • Среда без географии объясняет 212.89
(I_env <- sum_prda_1$constr.chi)
## [1] 213

Решение 3) Инерция, объясненная географией

sum_prda_2 <- summary(bf_prda_2)
> sum_prda_2
partit(sum_prda_2)
## Partitioning of variance:
##               Inertia Proportion
## Total             730      1.000
## Conditioned       279      0.383
## Constrained       193      0.265
## Unconstrained     257      0.352

География без среды объясняет 193.24

(I_geo <- sum_prda_2$constr.chi)
## [1] 193

Решение: 4) Инерция, совместно объясненная средой и географией

(I_env_geo <- I_total - I_env - I_geo)
## [1] 66.5

Компоненты инерции — сводим результаты вместе

comp <- data.frame(Inertia = c(I_env, I_geo, I_env_geo, I_total))
rownames(comp) <- c('Только среда', 
                    'Только география', 
                    'Среда и география вместе', 
                    'Общая объяснимая инерция')
comp$Proportion <- comp$Inertia/sum(comp$Inertia[1:3]) * 100
colnames(comp) <- c('Инерция', '%')
comp
##                          Инерция     %
## Только среда               212.9  45.0
## Только география           193.2  40.9
## Среда и география вместе    66.5  14.1
## Общая объяснимая инерция   472.6 100.0

Среда объясняет 45% общей изменчивости генетической структуры — очень много, но и география объясняет 41%. И только 14% объясняется совместным влиянием среды и географии

Take home messages

  • Анализ избыточности помогает установить связь между несколькими наборами переменных. Один из наборов считается зависимым, другой считается объясняющим
  • Для анализа необходимо, чтобы зависимости переменных-откликов от предикторов были линейными
  • В ходе анализа выделяют два типа осей — канонические (ограниченные, объясненные) переменными-предикторами, и неканонические (неограниченные, необъясненные) ими
  • Частный анализ избыточности позволяет описать зависимость двух наборов переменных с поправкой на влияние дополнительных переменных (ковариат)
  • При помощи частного анализа избыточности можно выделить компоненты инерции связанные с несколькими (2–4) наборами переменных-предикторов

Дополнительные ресурсы

  • Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.
  • Legendre, P., Legendre, L., 2012. Numerical ecology. Elsevier.
  • Oksanen, J., 2011. Multivariate analysis of ecological communities in R: vegan tutorial. R package version 2–0.
  • The Ordination Web Page URL http://ordination.okstate.edu/ (accessed 19.04.17).
  • Quinn, G.G.P., Keough, M.J., 2002. Experimental design and data analysis for biologists. Cambridge University Press.