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

library(vegan)
data(mite)
data(mite.env)
data(mite.xy)
head(mite[ , 1:6], 2)

##   Brachy PHTH HPAV RARD SSTR Protopl
## 1     17    5    5    3    2       1
## 2      2    7   16    0    6       0

str(mite.xy)

## 'data.frame':    70 obs. of  2 variables:
##  $ x: num  0.2 1 1.2 1.4 2.4 1.8 0.05 2 2 1.2 ...
##  $ y: num  0.1 0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5 1.7 ...

str(mite.env)

## 'data.frame':    70 obs. of  5 variables:
##  $ SubsDens : num  39.2 55 46.1 48.2 23.6 ...
##  $ WatrCont : num  350 435 372 360 204 ...
##  $ Substrate: Factor w/ 7 levels "Sphagn1","Sphagn2",..: 1 5 7 1 1 1 1 7 5 1 ...
##  $ Shrub    : Ord.factor w/ 3 levels "None"<"Few"<"Many": 2 2 2 2 2 2 2 3 3 3 ...
##  $ Topo     : Factor w/ 2 levels "Blanket","Hummock": 2 2 2 2 2 2 2 1 1 2 ...

library(ggplot2)
library(gridExtra)
th <- theme(legend.text = element_text(size = 13), legend.title = element_text(size = 14), legend.position = "bottom", text = element_text(size = 10))
th_narrow <- th + theme(legend.key.width = unit(3, "mm"))

theme_set(theme_bw(base_size = 18))
update_geom_defaults("point", list(shape = 19))
p <- ggplot(mite.xy, aes(x = x, y = y)) + geom_point()
p + coord_fixed() + th

• ищем ось макс. варьирования (анализ главных компонент (PCA), корреспондентный анализ (CA))
  
• затем строим регрессионную модель, описывающую связь PC или CA от интересующих нас предикторов

• Опасности:
• теряем связи переменных среды и отброшенных компонент высоких порядков.
• Выход использовать прямой градиентный анализ.
• PCA нужны линейные зависимости (т.к. матрица ковариаций или корреляций). Если связи нелинейны - искривление градиентов.
• Выход использовать другие расстояния (CA использует хи-квадрат).

Прямой градиентный анализ еще называют каноническим анализом.

Происхождение термина канонический

Канонической формой функции или выражения в математике называют такую упрощенную форму, к которой можно свести функцию или выражение без потери самой важной информации.

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

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

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

НО! Оси ординации определяются еще и переменными среды.

Переменные среды ограничивают оси, поэтому канонические оси еще называют ограниченными осями (constrained axis)

\[ \hat{y_i} = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \dots + \beta_px_{ip} \]

Здесь \(\hat{y}\) - предсказанные значения (fitted values)

\(x_1, x_2 \dots x_p\) - предикторы

\(\beta_0, \beta_1, \beta_2 \dots + \beta_p\) - коэффициенты

Как всегда, независимая часть модели - это модельная матрица X
Зависимая часть в анализе - это не вектор y, как мы привыкли, а матрица Y
Коэффициенты, будут представлены тоже в виде матрицы \(\beta\).

Тогда матрица предсказанных значений

\[ \hat {\textbf Y} = \boldsymbol {\beta \times X} \]

При этом матрица коэффициентов - это

\[ \boldsymbol{\beta} = \boldsymbol {[X'X]^{-1}X'Y} \]

Матрица остатков - это \[ \boldsymbol {Y_{resid} = Y - \hat {Y}} \]

Таким образом на вход в канонический анализ поступают три матрицы: \(\boldsymbol{Y}\), \(\boldsymbol{\hat{Y}}\), \(\boldsymbol{Y_{resid}}\)

1. SVD матрицы исходных значений (\(\boldsymbol{Y}\)) - дает информацию об общей изменчивости в системе - это просто СА (или PCA).
   
2. SVD матрицы предсказанных значений (\(\boldsymbol{\hat{Y}}\)) - дает информацию о системе, ограниченной (Constrained) предикторами (\(\boldsymbol{X}\)).
   
3. SVD матрицы остатков (\(\boldsymbol{Y_{resid}}\)) - дает информацию о системе за вычетом влияния предикторов.

• Основан на корреспондентном анализе, то есть вместо исходных сырых, центрированных данных (\(\boldsymbol{Y}\)) используется матрица \(\boldsymbol{Q}\) (Матрица вкладов в отклонение от нулевой модели, основанной на независимости между признаками и объектами).

Условия применимости такие же как у CA

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

Рассмотрение ординации в канонических осях дает информацию о связи между видами, сайтами и средой.

Анализ неканонических осей - отношения между видами или сайтами после того, как удален эффект среды.

• Зависимые переменные (отклики) - обилие видов

• Независимые переменные (предикторы) - переменные среды

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

В качестве предикторов в модели возьмите следующие факторы из датафрейма mite.env: SubsDens, WatrCont, Substrate, Topo

mite_cca <- cca(mite ~ SubsDens + WatrCont + Substrate + Topo, data = mite.env)

vif.cca(mite_cca)

##           SubsDens           WatrCont   SubstrateSphagn2 
##               1.43               1.93               1.16 
##   SubstrateSphagn3   SubstrateSphagn4    SubstrateLitter 
##               1.29               1.06               1.22 
##  SubstrateBarepeat SubstrateInterface        TopoHummock 
##               1.05               1.61               1.72

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

Partitioning of mean squared contingency coefficient:
              Inertia Proportion
Total           1.696      1.000
Constrained     0.735      0.433
Unconstrained   0.961      0.567

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

Accumulated constrained eigenvalues
Importance of components:
                       CCA1  CCA2   CCA3   CCA4   CCA5   CCA6 ...
Eigenvalue            0.426 0.129 0.0667 0.0443 0.0321 0.0145 ...
Proportion Explained  0.580 0.175 0.0907 0.0603 0.0437 0.0197 ...
Cumulative Proportion 0.580 0.755 0.8458 0.9061 0.9498 0.9695 ...

Сколько всего осей?

• Почему неканонических осей 34, а не 35?

• А почему канонических осей 9?

\[ \textbf{Q} = \frac {p_{ij} - p_ip_j}{\sqrt{p_ip_j}} = \frac{f_{ij}Ft - f_if_j}{Ft \sqrt{f_if_j}} \]

#Матрица вкладов, вычисленная через частоты
Q <- (f_ij*Ft - f_i %*% t(f_j))/(Ft*sqrt(f_i %*% t(f_j))) 
Q <- as.matrix(Q)

\[ \textbf{Q}_{n \times p} = \textbf{U}_{n \times p} \textbf{D}_{p \times p} \textbf{V}'_{p \times p} \]

# SVD матрицы Q
U <- svd(Q)$u 
D <- diag(svd(Q)$d)
V <- svd(Q)$v

Inertia_total <- sum(D^2) #Общая инерция в системе
CA_number <- sum(round(D, 2) !=0) #Количество главныю осей в CA. Обратите внимание на то, что их на одну меньше, чем исходных колонок в матрице Q. 

Эта величина оценивается “инерцией”

Общая изменчивость в системе 1.696

Общее количество неканонических (unconstrained) осей 34, как и положено, на одно меньше, чем общее количество колонок в матрице Q, равное 35

НО! Почему ограниченных осей 9, хотя предикторов у нас всего 4?

Для этого надо будет решить следующее матричное уравнение

\[ \hat {\textbf Y} = \boldsymbol {\beta \times X} \]

В нашем случае вместо \(\hat {\textbf Y}\) будет матрица Q_pred

\[ \boldsymbol{\beta} = \boldsymbol {[X'X]^{-1}X'Y} \]

\[ \boldsymbol {Y_{resid} = Y - \hat {Y}} \]

В нашем случае это будет матрица Q_resid

Поскольку строки и столбцы матрицы Q равнозначны, то есть два равнозначных центроида: центроид строк (пробы) и центроид колонок(виды)

Координаты центроидов задаются векторами p_i(для строк) и p_j(для колонок)

Для работы со строками надо учитывать вектор p_i

Для работы со столбцами надо учитывать вектор p_j

Матрица коэффициентов уравнения регрессии

\[ \boldsymbol{\beta} = \boldsymbol {[X'X]^{-1}X'Y} \\ \boldsymbol{\beta} = \boldsymbol{[X' p_i X]^{-1}[X'\sqrt{p_i} Q]} \]

#Матрица коэффициентов
betas <- solve(t(X) %*% diag(p_i) %*% X) %*% (t(X) %*% diag(p_i)^(1/2) %*% Q)

Обратите внимание, что для вычисления коэффициентов регрессии нам пришлось опять вводить в расчет маргинальные вероятности \(p_i\).

Это связано с природой матрицы Q. Эта матрица не является матрицей первичных данных, а является матрицей остатков от нулевой модели, описывающей отсутствие связи между столбцами и сроками. Поэтому все ее значения имеют смысл только по отношению к центроиду строк или по отношению к центроиду колонок.

Матрица предсказанных значений

\[ \boldsymbol{Q_{pred}} = \boldsymbol{\sqrt{p_i} X \beta} \]

 #Матрица предсказанных значенй
Q_pred <- diag(p_i)^(1/2) %*% X %*% betas

SVD матрицы предсказанных значений

U_pred <- svd(Q_pred)$u 
D_pred <- diag(svd(Q_pred)$d)
V_pred <- svd(Q_pred)$v

D_pred - характеризует канонические (constrained) оси

Inertia_constrained <- sum(D_pred^2) #Инерция в ограниченной ординации

CCA_number <- sum(round(D_pred, 2) !=0) #Количество Канонических осей в СCA. 

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

Инерция в канонических осях составляет 0.735.

Можно оценить, сколько процентов составляет инерция канонических осей от общей инерции, полученной в CA (1.696).

Итого, канонические оси отражают 43.343 % общей изменчивости (инерции).

Это разность между матрицей Q и Q_pred

Q_resid <- Q - Q_pred

U_res <- svd(Q_resid)$u 
D_res <- diag(svd(Q_resid)$d)
V_res <- svd(Q_resid)$v

Сингулярные числа в матрице D_res задают неканонические (unconstrained) оси

Число ненулевых сингулярных чисел в матрице D_res составляет 34. То есть неканонических осей столько же сколько главных осей в корреспондентном анализе. То есть, на 1 меньше, чем в колонок в матрице Q.

#Инерция ординации остатков 
Inertia_unconstrained <- sum(D_res^2) 

Инерция в пространстве неканонических осей составляет: 0.961, то есть неканонические оси объясняют 56.657 % общего варьирования (общей инерции) в системе.

На всякий случай обратите внимание, что \[Inertia_{total} = Inertia_{constrained} + Inertia_{unconstrained}\]

Inertia_unconstrained + Inertia_constrained

## [1] 1.7

Вот то, что показывает cca()

mite_cca <- cca(mite ~ SubsDens + WatrCont + Substrate + Topo, data = mite.env)
mite_cca

## Call: cca(formula = mite ~ SubsDens + WatrCont + Substrate +
## Topo, data = mite.env)
## 
##               Inertia Proportion Rank
## Total           1.696      1.000     
## Constrained     0.735      0.433    9
## Unconstrained   0.961      0.567   34
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##  CCA1  CCA2  CCA3  CCA4  CCA5  CCA6  CCA7  CCA8  CCA9 
## 0.426 0.129 0.067 0.044 0.032 0.014 0.011 0.008 0.004 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.1358 0.1205 0.1127 0.0916 0.0639 0.0575 0.0451 0.0404 
## (Showing 8 of 34 unconstrained eigenvalues)

#Координаты видов канонических осях полученные вручную

constr_CA_species <- diag(p_j^(-1/2))%*% V_pred

Виды – это колонки, стало быть, маргинальная вероятность p_j

#Координаты видов в канонических осях по версии cca()

cca_sp_constr <- mite_cca$CCA$v 

• Часть изменчивости структуры сообществ орибатид действительно удается объяснить изменениями среды. Канонические оси объясняют 43% общей инерции если рассматривать всю изменчивость (т.е. 43 оси). Из них первые две оси объясняют 33% общей изменчивости.

Но, если рассматривать первые две канонические оси в пределах только ограниченных осей (т.е. 9 осей) - то учтенные факторы среды объясняют целых 75%

Accumulated constrained eigenvalues
Importance of components:
                       CCA1  CCA2   CCA3   CCA4   CCA5   CCA6   CCA7 ...
Eigenvalue            0.426 0.129 0.0667 0.0443 0.0321 0.0145 0.0107 ...
Proportion Explained  0.580 0.175 0.0907 0.0603 0.0437 0.0197 0.0145 ...
Cumulative Proportion 0.580 0.755 0.8458 0.9061 0.9498 0.9695 0.9840 ...

• Первая ограниченная ось объясняет 58%, вторая 17% того, что можно потенциально связать с воздействием среды, остальные оси почти ничего не объясняют. Значит, характеристики среды, в принципе, можно свести к двум комплексным независимым переменным.

• обилие каких видов сильнее варьирует вдоль оси?

scores(mite_cca, display = "species", choices = 1:5)

• Собственные векторы (“species scores”) координаты видов
• Координаты объектов (“site scores”)
• Координаты объектов в канонических осях (“Site constraints”)
• (“Biplot scores”) координаты для биплотов
• Центроиды сайтов для бинарных переменных среды на диаграмме ординации (“Centroids for factor constraints”)

spenvcor(mite_cca)

##  CCA1  CCA2  CCA3  CCA4  CCA5  CCA6  CCA7  CCA8  CCA9 
## 0.905 0.771 0.641 0.697 0.605 0.527 0.521 0.448 0.269

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

Триплоты:

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

Биплоты:

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

Осторожно, в CCA есть два типа координат объектов:

• WA - “weighted average scores” - без ограничений переменными среды, но при этом они отличаются от координат полученных в CA

• LC - “linear combination scores” - результат множественной регрессии WA-координат по линейным комбинациям переменных среды - это по-умолчанию используется в vegan

• Тестируем гипотезу о том, что отношения между структурой сообщества и средой значимы - \(H _0\): обилие видов в пробах не зависит от значений переменных среды

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

• статистика - сумма всех канонических собственных чисел

Есть ли связь структуры сообщества со средой?

• Может основываться только на неотрицательных величинах.
• Редкие виды - использовать когда репрезентативная выборка или удалить из анализа
• Чтобы снизить вес больших численностей видов - логарифмировать
• Доля объясненной изменчивости (% of total inertia) - аналогично \(R^2\) - смещенная оценка. Простых поправок не придумано.
• Симметричные распределения переменных среды (преобразования)
• Нельзя включать шумные переменные среды

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

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

• Вычислительная часть CCA ничем не отличается вычислительной части CA, но в анализ вовлекаются вместо одной матрицы Q две матрицы:

• Матрица предсказанных линейной моделью значений Q_pred (дает канонические оси)

• Матрица остатков Q_res (дает неканонические оси)

• Вычисление матриц Q_pred и Q_res производится с использованием матричного способа подбора коэффициентов регрессии.

• Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.
• Jongman, R.H.G., Ter Braak, C.J.F., Van Tongeren, O.F.R. (eds.), 1995. Data analysis in community and landscape ecology. Cambridge University Press
• 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 10.21.13).
• Quinn, G.G.P., Keough, M.J., 2002. Experimental design and data analysis for biologists. Cambridge University Press.