class: middle, left, inverse, title-slide .title[ # RDA: redundancy analysis, анализ избыточности ] .subtitle[ ## Анализ и визуализация многомерных данных с использованием R ] .author[ ### Марина Варфоломеева ] .author[ ### Вадим Хайтов ] .author[ ### Анастасия Лянгузова ] --- ## Анализ избыточности (Redundancy analysis, RDA) - Связь нескольких наборов переменных - Анализ избыточности, теория и практика - Проверка значимости ординации - Выбор оптимальной модели - Частный анализ избыточности и компоненты объясненной инерции - Компоненты объясненной изменчивости ### Вы сможете - Проводить анализ избыточности - Оценивать долю объясненной инерции - Интерпретировать компоненты по нагрузкам переменных - Строить ординацию объектов в пространстве компонент - Проверять значимость модели ординации при помощи пермутационного теста - Разделять объясненную инерцию на компоненты, связанные с разными наборами переменных, при помощи частного анализа избыточности --- class: middle, center, inverse # Связь нескольких наборов переменных --- ## Что будет определять генетическую структуру в колониях бабочек? ### Пример: генетика бабочек _Euphydryas editha_ Частоты разных аллелей фосфоглюкоизомеразы и данные о факторах среды для 16 колоний бабочек _Euphydryas editha_ в Калифорнии и Орегоне (данные McKechnie et al., 1975) .pull-left[. <!-- --> ] .pull-right[  ] <!-- фосфоглюкоизомераза превращает глюкозу во фруктозу, участвует в гликолизе и в глюконеогенезе. нейротрофический фактор для спинальных и чувствительных нейронов, и т.п. --> --- class: middle, center, inverse # Анализ избыточности --- ## Анализ избыточности (Redundancy analysis, RDA) .content-box-green[ RDA --- метод прямой ординации из группы методов канонического анализа (ограниченной ординации = constrained ordination). ] RDA совмещает в себе: - множественную линейную регрессию - анализ главных компонент (PCA, principal component analysis) Особенности: - две матрицы данных: матрица предикторов (X, размерность n × m) и матрица зависимых переменных/переменных отклика (Y, размерность n × p); - поиск компонент из матрицы зависимых переменных, которые являются линейными комбинациями предикторов и отражают максимум изменчивости. --- ## Ординации и регрессия <img src="images/rda_ordination.png" height=300px> -- ![:col_row PCA --- анализ главных компонент, Множественная регрессия, RDA --- анализ избыточности] ![:col_row CA --- корреспондентный анализ, , CCA --- канонический корреспондентный анализ] Рис. из Legendre, Legendre, 2012 с изменениями --- ## 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]$$` --- ## Математика RDA 1. Множественная линейная регрессия для каждой зависимой переменной. Расчёт коэффициентов: `$$b = [\mathbf{X}'\mathbf{X}]^{-1} \mathbf{X}'\mathbf{y}$$` Можно получить матрицу __всех__ коэффициентов регрессии: `$$B = [\mathbf{X}'\mathbf{X}]^{–1} \mathbf{X}'\mathbf{Y}$$` -- Размерность такой матрицы будет m x p. --- ## Математика RDA 2. Расчёт предсказанных значений. Во множественной линейной регрессии считается как: `$$\hat{y} = \mathbf{X}\mathbf{b} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}(\mathbf{X}'\mathbf{y})$$` -- В случае RDA зависимая переменная --- матрица, которая будет считаться как: `$$\mathbf{\hat Y} = \mathbf{X}\mathbf{b} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}(\mathbf{X}'\mathbf{Y})$$` --- ## Условия применимости RDA Сходны с линейной регрессией. -- 1. Независимость наблюдений. 2. Линейная зависимость переменных-откликов от предикторов. 3. Отсутствие мультиколлинеарности предикторов. 4. Наблюдений должно быть __значительно больше__, чем предикторов (как в регрессии, см. "проклятие размерности") --- ## Подготовка данных к RDA 1. Центрирование матриц `\(\mathbf{Y}\)` и `\(\mathbf{X}\)`. При необходимости --- стандартизация. 2. Проверка на мультиколлинеарность (как в линейной регрессии). 3. Проверка на нормальное распределение (как в линейной регрессии). --- ## Вычисление канонических осей .pull-left[ <div class="DiagrammeR html-widget html-fill-item" id="htmlwidget-7e94875635a0c8153be6" style="width:432px;height:540px;"></div> <script type="application/json" data-for="htmlwidget-7e94875635a0c8153be6">{"x":{"diagram":"graph TD; \n A(Центрирование матрицы Y)-->B(Множественная линейная регрессия)\n B-->C(Матрица предсказанных значений); \n C-->D(Анализ главных компонент);\n D-->E(Матрица собственных векторов); E-->F(Ординация в пространстве Y);\n E-->G(Ординация в пространстве X); A-->F; C-->G\n style A fill:#C6FFE7,stroke:#C6FFE7;\n style B fill:#C6FFE7,stroke:#C6FFE7;\n style C fill:#C6FFE7,stroke:#C6FFE7;\n style D fill:#C6FFE7,stroke:#C6FFE7;\n style E fill:#C6FFE7,stroke:#C6FFE7;\n style F fill:#C6FFE7,stroke:#C6FFE7;\n style G fill:#C6FFE7,stroke:#C6FFE7;\n "},"evals":[],"jsHooks":[]}</script> ] .pull-right[1. Множественная линейная регрессия `\(\mathbf{Y}_{np}\)` от `\(\mathbf{X}_{nm}\)`. Получаем матрицу предсказанных значений: `\(\mathbf{\hat Y} = \mathbf{X}[\mathbf{X}'\mathbf{X}]^{-1}\mathbf{X}'\mathbf{Y}\)` 2. PCA по `\(\mathbf{\hat Y}\)`. Получаем собственные числа и векторы (матрица `\(\mathbf{U}\)`). - Ординация объектов в пространстве наблюдаемых значений `\(\mathbf{Y}\)`: `$$\mathbf{F} = \mathbf{Y}\mathbf{U}$$` - Ординация объектов в пространстве предикторов `\(\mathbf{X}\)`: `$$\mathbf{Z} = \mathbf{\hat{Y}} \mathbf{U}$$`] --- ## Вычисление неканонических осей .pull-left[ <div class="DiagrammeR html-widget html-fill-item" id="htmlwidget-ebf03018d20da17db214" style="width:432px;height:540px;"></div> <script type="application/json" data-for="htmlwidget-ebf03018d20da17db214">{"x":{"diagram":"graph TD; \n A(Получение остатков от множественных регрессий)-->B(Анализ главных компонент)\n B-->C(Матрица собственных векторов); \n C-->D(Ординация в пространстве остатков);\n A-->D\n style A fill:#C6FFE7,stroke:#C6FFE7;\n style B fill:#C6FFE7,stroke:#C6FFE7;\n style C fill:#C6FFE7,stroke:#C6FFE7;\n style D fill:#C6FFE7,stroke:#C6FFE7;\n "},"evals":[],"jsHooks":[]}</script> ] .pull-right[ 3. PCA по остаткам от регрессии: `\(\mathbf{Y_{res}} = \mathbf{Y - \hat Y}\)`. Получаем собственные векторы остатков `\(\mathbf{U_{res}}\)`. - Ординация объектов в пространстве остатков от регрессии: `$$\mathbf{Y_{res}U_{res}}$$` ] --- class: middle, center, inverse # Подготовка данных к RDA --- ## Рассмотрение исходных данных Нам нужны две матрицы данных: матрица предикторов и матрица зависимых переменных/переменных отклика. Структура исходных данных по 16 колониям. - `$xy` --- координаты колоний - `$envir` --- 4 фактора среды для колоний - `$genet` --- частоты 6 аллелей в колониях - `$contour` --- карта Калифорнии - `$Spatial` --- объект, содержащий карту. ```r library(ade4) data(butterfly) str(butterfly) ``` --- ## Создадим переменные с более короткими названиями для удобства ```r gen <- butterfly$genet head(gen, 1) ``` ``` 0.4 0.6 0.8 1 1.16 1.3 SS 0 3 22 57 17 1 ``` ```r 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 ``` --- class: middle, center, inverse # Проверка условий применимости --- ## Стандартизация данных ```r summary(gen) ``` ``` 0.4 0.6 0.8 1 1.16 Min. : 0.00 Min. : 0.00 Min. : 1.0 Min. :25.0 Min. : 0.0 1st Qu.: 0.00 1st Qu.: 2.75 1st Qu.:12.5 1st Qu.:36.8 1st Qu.:10.0 Median : 0.00 Median : 4.50 Median :18.0 Median :47.0 Median :17.0 Mean : 1.75 Mean : 7.19 Mean :18.6 Mean :51.2 Mean :17.2 3rd Qu.: 0.25 3rd Qu.: 7.50 3rd Qu.:23.5 3rd Qu.:59.2 3rd Qu.:27.0 Max. :14.00 Max. :26.00 Max. :40.0 Max. :92.0 Max. :35.0 1.3 Min. : 0.00 1st Qu.: 0.00 Median : 3.00 Mean : 4.12 3rd Qu.: 6.50 Max. :14.00 ``` Стандартизация не нужна. --- ## Линейность связи и мультиколлинеарность ```r pairs(env_geo) ``` <!-- --> --- ## Удаляем коллинеарные предикторы Считаем, что VIF `\(\geq\)` 10 --- значительная коллинеарность (Borcard et al., 2011). ```r library(car) vif(lm(gen$`0.4` ~ ., data = env_geo)) ``` ``` Altitude Precipitation Temp_Max Temp_Min x y 23.889 4.456 4.769 22.492 7.544 11.245 ``` -- - Минимальная температура связана с высотой. Придется оставить что-то одно. Не будем использовать минимальную температуру. ```r vif(lm(gen$`0.4` ~ . -Temp_Min, data = env_geo)) ``` ``` Altitude Precipitation Temp_Max x y 7.516 4.407 4.606 5.612 7.271 ``` --- class: middle, center, inverse # RDA в R --- ## RDA в vegan - Зависимые переменные (отклики) --- генетические данные - Независимые переменные (предикторы) --- переменные среды ```r 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 PC4 PC5 Eigenvalue 352.481 18.4680 0.78079 212.437 108.097 27.2173 8.9505 1.1684 Proportion Explained 0.483 0.0253 0.00107 0.291 0.148 0.0373 0.0123 0.0016 Cumulative Proportion 0.483 0.5084 0.50950 0.801 0.949 0.9861 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.23 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 729.6 1.000 Constrained 371.7 0.509 Unconstrained 357.9 0.491 ``` Total --- для всех осей, общая изменчивость исходной матрицы откликов (генетической структуры в разных сайтах); Constrained --- для осей, которые являются комбинациями предикторов (изменчивость объясненная средой); Unconstrained --- необъясненная предикторами изменчивость. --- ## Влияние компонент Можно более подробно оценить, как распределяется изменчивость между осями ``` Eigenvalues, and their contribution to the variance ``` ``` $importance Importance of components: RDA1 RDA2 RDA3 PC1 PC2 PC3 PC4 PC5 Eigenvalue 352.481 18.4680 0.78079 212.437 108.097 27.2173 8.9505 1.1684 Proportion Explained 0.483 0.0253 0.00107 0.291 0.148 0.0373 0.0123 0.0016 Cumulative Proportion 0.483 0.5084 0.50950 0.801 0.949 0.9861 0.9984 1.0000 ``` Больше всего объясняет первая каноническая ось (48.3%), в то же время первые две неканонические --- ещё 43.9%. --- ## Распределение изменчивости, потенциально объяснимой факторами ``` 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” ```r scores(bf_rda, display = "species", choices = 1:5) ``` ``` RDA1 RDA2 RDA3 PC1 PC2 0.4 0.6786 -0.57904 0.21624 -0.9465 -0.55116 0.6 0.9097 -0.73038 -0.12097 -2.2304 -0.76028 0.8 2.6171 -0.05586 -0.18522 -2.7724 0.03138 1 -6.3060 0.11997 -0.05311 3.3960 -2.50904 1.16 1.5146 1.32492 0.03152 2.3040 2.65044 1.3 0.5860 -0.07961 0.11153 0.2492 1.13867 attr(,"const") [1] 10.23 ``` --- ## Корреляции между откликами и предикторами Сильная корреляция между генетической структурой и средой только для первой ограниченной оси. Для других --- умеренные или слабые. ```r spenvcor(bf_rda) ``` ``` RDA1 RDA2 RDA3 0.8304 0.3480 0.1663 ``` Будьте осторожны с интерпретацией! Возможна ситуация, когда корреляция между откликами и предикторами будет сильной, а соответствующие канонические оси будут объяснять ничтожную долю изменчивости. --- class: middle, center, inverse # Визуализация ординации --- ## Визуализация ординации - Какие предикторы важнее всего? - Какими факторами определяется значение зависимых переменных? ### Биплоты: - отклики + предикторы - объекты + предикторы ### Триплоты: - переменные-отклики ("species"), - объекты ("sites") - переменные-предикторы (непрерывные в виде векторов, дискретные в виде центроидов) --- ## Триплот корреляций (scaling = 2): Какие переменные среды сильнее всего определяют сходство объектов? .pull-left[ ```r plot(bf_rda, scaling = 2) ``` <!-- --> ] .pull-right[ - Векторы --- независимые переменные, факторы среды - Надписи --- объекты (сайты, особи, популяции и пр.) - Красные надписи --- зависимые переменные - Косинусы углов между векторами --- корреляции между соотв. переменными - Расстояния между объектами не имеют конкретного смысла - Проекция объекта на линию-вектор --- значение переменной для данного объекта ] --- ## Пример интерпретации триплота корреляций ```r plot(bf_rda, scaling = 2) ``` <!-- --> - Вдоль первой оси изменяется температура, высота и осадки - Вдоль второй оси --- немного меняется уровень осадков --- ## Триплот расстояний (scaling = 1): Насколько похожи друг на друга объекты? .pull-left-33[ ```r plot(bf_rda, scaling = 1) ``` <!-- --> ] .pull-right-66[ Надписи --- объекты (сайты, особи, популяции и пр.). Красные надписи --- зависимые переменные. Векторы --- независимые переменные, факторы среды. Расстояния между точками --- расстояния между наблюдениями. Косинусы углов между векторами предикторов и откликов --- корреляции между соотв. переменными, другие углы не имеют смысла. Проекция объекта на линию-вектор отражает примерное относительное положение данного объекта вдоль соответствующей переменной (но не значение). Отношения между дискретными и непрерывными предикторами не интерпретируются. ] --- ## Пример интерпретации триплота расстояний .pull-left[ ```r plot(bf_rda, scaling = 1) ``` <!-- --> ] .pull-right[ - Генетическая структура в LO и UO похожа, но не похожа на остальные места - GL и GH --- более высокогорные сайты, чем LO и UO ] --- class: middle, center, inverse # Проверка значимости ординации --- ## Оценка значимости полученной ординации 1. Существует ли зависимость от предикторов. 2. Влияние факторов на зависимые переменные. 3. Поиск конкретных факторов, влияющих на переменные (эффекты type I и III). 4. Значимость изменчивости вдоль ординационных осей. --- ## R-статистика `\(\mathbf{R^2}\)` измеряет силу связи между `\(\mathbf{Y}\)` и `\(\mathbf{X}\)`. `\(\mathbf{R^2_{adj}}\)` делает то же самое, но учитывает число переменных в матрице предикторов `\(\mathbf{X}\)`. `$$\mathbf{R^2_{adj}} = 1 - (1 - \frac{\mathbf{SS(\hat{Y})}}{\mathbf{SS (Y)}}) \frac{n - 1}{n - m - 1}$$` `\(n\)` --- число наблюдений, `\(m\)` --- ранг канонических осей (число предикторов). --- ## Псевдо-F статистика Псевдо-F статистика расчитывает значимость полученной ординации. `$$\mathbf{pseudoF} = \frac{\mathbf{R^2_{adj}} \mathbf{m}}{\mathbf{(1 - R^2_{adj})(n-m-1)}}$$` `\(n\)` --- число наблюдений, `\(m\)` --- ранг канонических осей (число предикторов). Это же уравнение может быть переписано с использованием значений инерции. -- `$$\mathbf{pseudoF = \frac {\Lambda_{c} / m} {\Lambda_{r}/(n - m - 1)}}$$` `\(\Lambda_{c}\)` --- инерция канонических осей; `\(\Lambda_{r}\)` --- остаточная инерция. Псевдо-F статистика не следует F-распределению как в ANOVA. Чтобы оценить распределение, нужны пермутации. --- ## Пермутации -- 1. Пермутация (перестановка) строк матрицы зависимых переменных `\(\mathbf{Y}\)`. 2. Расчёт псевдо-F статистики для пермутации. 3. Повторяем пункты 1-2 много-много раз. По распределению пермутационной статистики оцениваем долю пермутаций, в которых качество модели оказалось лучше, чем для реальной модели --- это доверительная вероятность `\(p\)`. Если в большинстве пермутаций получаются модели хуже, чем реальная модель, то значит зависимость от предикторов значима. --- ## Общий тест: Влияют ли факторы на зависимые переменные? - тестируем гипотезу о том, что отношения между генотипом и средой значимы. `\(H _0\)`: значения откликов в пробах не зависят от переменных среды (генетическая структура не зависит от среды) ```r 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.0067 ** Residual 12 358 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Есть ли связь генетической структуры со средой? -- Связь генетической структуры и среды значима. --- ## Тест факторов, type I эффекты: Какие факторы влияют на зависимые переменные? ```r 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.0016 ** Precipitation 1 72 2.43 0.1052 Temp_Max 1 20 0.67 0.5109 Residual 12 358 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` -- Генетическая структура популяций бабочек достоверно зависит от высоты, если в модель включены др. факторы. Но это Type I эффекты --- они зависят от порядка включения факторов в модель. Т.е. после включения высоты в модель другие факторы уже не влияют. --- ## Тест факторов, type III эффекты: Какие факторы влияют на зависимые переменные? ```r 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.52 Residual 12 358 ``` -- Если протестировать каждый из факторов отдельно, при условии, что все остальные включены в модель, то получится, что ни один из них не влияет. --- ## Тест значимости осей, ограниченных факторами: Вдоль какой из осей значимо меняется генетическая структура? - `\(H _0\)`: значения переменных-откликов для объектов не зависят от переменных-предикторов - пермутационный: выбирает оси, которые объясняют больше изменчивости, чем из др. матриц, полученных путем перестановок ```r 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.8495 RDA3 1 1 0.03 0.9983 Residual 12 358 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` -- Генетическая структура значимо меняется вдоль первой главной оси. --- class: middle, center, inverse # Выбор оптимальной модели --- ## Выбор оптимальной модели В нашей модели много факторов. Если мы тестируем любой из факторов, после включения остальных в модель --- он не влияет. Вероятно, модель не оптимальна. Как подобрать оптимальную модель? -- > - Можно использовать пошаговый выбор модели: добавляем в модель лучшие переменные и снова исключаем те, что потеряли значимость. (Вспомните, как это было для регрессионных моделей.) Какой можно использовать тест для сравнения моделей? Модели с разным числом предикторов можно сравнить при помощи пермутационного теста (AIC для канонических ординаций не существует!) __Осторожно!__ Пошаговый выбор модели --- не панацея, т.к. разные пошаговые методы могут дать разные конечные модели! Не запихивайте в модель сразу все предикторы, а начинайте ваш выбор с небольшого числа важных факторов. Финальная модель не обязательно верна! --- ## Пошаговый выбор оптимальной модели Пошаговый выбор модели может быть: 1. "Backward selection" --- начинает работать с __полной__ модели, а затем постепенно удаляет наименее значимые один за другим. 2. "Forward selection" --- начинать работать с __нулевой__ модели, постепенно добавляя переменные, вносящие наибольший вклад. --- ## Пошаговый выбор оптимальной модели Для пошагового выбора нам понадобятся полная и нулевая модели ```r m1 <- rda(gen ~ Altitude + Precipitation + Temp_Max, data = env_geo) m0 <- rda(gen ~ 1, data = env_geo) ``` Запускаем пошаговый выбор. В случае указания scope направление пошаговой модели по умолчанию "both". ```r m <- ordistep(m0, scope = formula(m1), permutations = 9999) ``` ``` Start: gen ~ 1 Df AIC F Pr(>F) + Altitude 1 101 8.69 0.0022 ** + Temp_Max 1 101 8.12 0.0035 ** + 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.0018 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Df AIC F Pr(>F) + Precipitation 1 99.9 2.49 0.095 . + Temp_Max 1 102.0 0.64 0.551 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Оптимальная модель, отобранная при помощи пошагового алгоритма ```r m$anova ``` ``` Df AIC F Pr(>F) + Altitude 1 101 8.69 0.0022 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Оптимальная модель содержит только один предиктор --- высоту. --- ## Пошаговый выбор оптимальной модели Модель можно выбрать на основании значения `\(R^2_{adj}\)` и p-value. ```r m_radj <- ordiR2step(m0, scope = formula(m1), permutations = 9999) ``` ``` Step: R2.adj= 0 Call: gen ~ 1 R2.adjusted <All variables> 0.3869 + Altitude 0.3389 + Temp_Max 0.3219 + Precipitation 0.2958 <none> 0.0000 Df AIC F Pr(>F) + Altitude 1 101 8.69 0.002 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Step: R2.adj= 0.3389 Call: gen ~ Altitude R2.adjusted + Precipitation 0.4026 <All variables> 0.3869 <none> 0.3389 + Temp_Max 0.3214 ``` --- ## Оптимальная модель Также значимым предиктором является исключительно высота. ```r m_radj$anova ``` ``` R2.adj Df AIC F Pr(>F) + Altitude 0.339 1 101 8.69 0.002 ** <All variables> 0.387 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- class: middle, center, inverse # Частный анализ избыточности --- ## Применение частного RDA Так же как в случае обычного RDA есть две матрицы данных: матрица предикторов `\(\mathbf{X}\)` и матрица зависимых переменных/переменных отклика `\(\mathbf{Y}\)`, к ним добавляется матрица __ковариат__ `\(\mathbf{W}\)`. В случае с бабочками мы обнаружили связь генотипов со средой, но эта связь может быть обусловлена разными причинами. -- ### Что делать? Нужно удалить влияние географического положения, чтобы сделать корректный вывод о связи генотипов со средой. Тут нам и пригодится матрица ковариат. --- ## Частный RDA Рассчитываем зависимость от одного набора переменных (предикторов), исключая влияние ковариат. ### Алгоритм <div class="DiagrammeR html-widget html-fill-item" id="htmlwidget-103dbda9ce16bfffd5f7" style="width:648px;height:360px;"></div> <script type="application/json" data-for="htmlwidget-103dbda9ce16bfffd5f7">{"x":{"diagram":"graph TD; \n A(Множественная регрессия зависимости предикторов от ковариат)-->B(Остатки Y по W и X по W)\n B-->C(RDA остатков Y и остатков X); \n B-->D(RDA исходной матрицы Y и остатков X);\n style A fill:#C6FFE7,stroke:#C6FFE7;\n style B fill:#C6FFE7,stroke:#C6FFE7;\n style C fill:#C6FFE7,stroke:#C6FFE7;\n style D fill:#C6FFE7,stroke:#C6FFE7;\n "},"evals":[],"jsHooks":[]}</script> --- ## Частный RDA в R Осуществляется функцией в пакете vegan. Примерный алгоритм: (из Legendre, Legendre, 2012). ```r pRDA <- function(Y, X = NULL, W = NULL, scale.Y = FALSE) { Y <- scale(as.matrix(Y), center = TRUE, scale = scale.Y) if (!is.null(W)) { # При наличии ковариат W W <- scale(as.matrix(W), center = TRUE, scale = FALSE) Y <- qr.resid(qr(W), Y) } if (!is.null(X)) { # При начилии матрицы X X <- scale(as.matrix(X), center = TRUE, scale = FALSE) X <- cbind(X, W) Q <- qr(X) RDA <- svd(qr.fitted(Q, Y)) RDA$w <- Y %*% RDA$v %*% diag(1/RDA$d) Y <- qr.resid(Q, Y) } else { # Если нет ни X, ни W RDA <- NULL } RES <- svd(Y) # Анализ главныз компонент по остаткам list(RDA = RDA, RES = RES) } ``` --- ## Частный RDA по бабочкам ```r 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.0013 ** Residual 12 257 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Высота объясняет генетическую изменчивость, даже после удаления влияния географической близости. --- ## График ординации ```r op <- par(mfrow = c(1, 2)) plot(bf_prda_1, main = "Correlation triplot", scaling = 2) plot(bf_prda_1, main = "Distance triplot", scaling = 1) ``` <!-- --> ```r par(op) ``` --- class: middle, center, inverse # Компоненты объясненной инерции --- ## Компоненты инерции Среда объясняет генетическую изменчивость, даже после удаления влияния географических координат. Но какая часть изменчивости генетической структуры объясняется в чистом виде географической близостью, а какая --- общим действием среды и географии? --- ## Общую инерцию делим на части .center[ <!-- --> ] | Компонент инерции | Что описывает? | Инерция | | ---- | ---- | ---- | | a + b + c | вся потенциально объяснимая средой и географией изменчивость | Объясненная каноническими осями RDA среда + география | | a | изменчивость, объясненная средой | Объясненная каноническими осями pRDA среда + Condition (география) | | b | изменчивость, объясненная географией | Объясненная каноническими осями pRDA география + Condition (среда) | | c | изменчивость, совместно объясненная средой и географией | Находим инерцию по разности | --- ## Подбираем модели RDA, нужные для поиска компонентов инерции .center[ <!-- --> ] Нам нужна __полная модель RDA: генетика от среды и географического положения__ (чтобы найти a + b + c) У нас уже есть __частный RDA №1: зависимость генетики от среды с учетом географии__ (чтобы найти a) Нам нужен __частный RDA №2: генетика от географии с учетом свойств среды__ (чтобы найти b) ```r 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) Сколько инерции потенциально объясняется средой и географией? ```r sum_full <- summary(bf_rda_full) ``` ``` Partitioning of variance: ``` ``` Inertia Proportion Total 729.6 1.000 Constrained 472.6 0.648 Unconstrained 257.0 0.352 ``` Инерция, объясненная вместе средой и географией, здесь достаточно велика --- 472.64 ```r (I_total <- sum_full$constr.chi) ``` ``` [1] 472.6 ``` В отличие от нее, доля инерции, объясненной ограниченной матрицей, может быть довольно малой по отношению к общей инерции, поэтому можно сосредоточиться на доле от потенциально объяснимой инерции (от `sum_full$constr.chi`) --- ## Решение 2. Инерция, объясненная средой ```r sum_prda_1 <- summary(bf_prda_1) ``` ``` Partitioning of variance: ``` ``` Inertia Proportion Total 729.6 1.000 Conditioned 259.8 0.356 Constrained 212.9 0.292 Unconstrained 257.0 0.352 ``` - Среда без географии объясняет 212.89 ```r (I_env <- sum_prda_1$constr.chi) ``` ``` [1] 212.9 ``` --- ## Решение 3. Инерция, объясненная географией ```r sum_prda_2 <- summary(bf_prda_2) ``` ``` Partitioning of variance: ``` ``` Inertia Proportion Total 729.6 1.000 Conditioned 279.4 0.383 Constrained 193.2 0.265 Unconstrained 257.0 0.352 ``` География без среды объясняет 193.24 ```r (I_geo <- sum_prda_2$constr.chi) ``` ``` [1] 193.2 ``` --- ## Решение 4. Инерция, совместно объясненная средой и географией ```r (I_env_geo <- I_total - I_env - I_geo) ``` ``` [1] 66.5 ``` --- ## Компоненты инерции --- сводим результаты вместе ```r 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.04 Только география 193.2 40.89 Среда и география вместе 66.5 14.07 Общая объяснимая инерция 472.6 100.00 ``` Среда объясняет 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.