• Сложности при анализе видового состава сообществ при помощи анализа главных компонент
  • Анализ сырых данных
• Корреспондентный анализ
  • Анализ таблиц сопряженности, хи-квадрат
  • Оси в корреспондентном анализе
  • Интерпретация графиков в корреспондентном анализе

Вы сможете

• Находить проявление “эффекта подковы” в анализе главных компонент
• Проводить корреспондентный анализ таблиц сопряженности
• Объяснить, что именно означает взаиморасположение точек объектов и переменных на графиках результатов корреспондентного анализа
• Интерпретировать графики результатов корреспондентного анализа

Обилие 102 видов птиц в 37 сайтах в юго-восточной Австралии (Mac Nally, 1989; данные из Quinn, Keough, 2002). Можно ли описать отношения между сайтами небольшим числом главных компонент?

library(readxl)
birds <- read_excel(path = "data/macnally.xlsx")

# имена переводим в нижний регистр
colnames(birds) <- tolower(colnames(birds))

Построим искусственный градиент

• Линейные связи между переменными (т.к. матрица корреляций или ковариаций, подразумевает линейность связей)
• Исключить наблюдения, в которых есть пропущенные значения
• Если много нулей - трансформация данных
• Если очень много нулей - удалить такие переменные из анализа (иногда это большая редукция!)
• Если есть один или несколько пространственных градиентов, то часто лучше отказаться от PCA и использовать другие типы анализов.

mite_pca <- rda(mite)

screeplot(mite_pca, bstick = T)

biplot(mite_pca,  scaling = "sites", type = "t")

mite_mds <- metaMDS(mite, trace = F)

Проблемы PCA:

• Есть отскакивающее значение
• Вырисовывается “подкова”

Pro

• MDS может работать с любыми данными
• MDS лучше отражает градиенты

Contra

• MDS не дает в явном виде информации о связи признаков
  
• Оси MDS безлики их информативность не поддается оценке
  
• Оси MDS трудно использовать как комплексные признаки

• Существует третий путь - Корреспондентный анализ

op <- par(mfrow = c(1, 3))
plot(mite_pca, display = "sites", main = "PCA", type = "p")
plot(mite_mds, display = "sites",  main = "MDS")
plot(mite_ca, display = "sites", type = "p", main = "CA")
par(op)

op <- par(mfrow = c(1, 2))

plot(mite_ca, display = "sites")
plot(mite_ca, display = "species")

par(op)

Можно на одном графике

op <- par(mfrow = c(1, 2))

plot(mite_ca, scaling = 2, main = "scaling = 1 (sites)")
plot(mite_ca, scaling = 1, main = "scaling = 2 (species)")
# plot(mite_ca, scaling = 1, main = "scaling = 3 ")

par(op)

В основе CA лежит исследование таблиц сопряженности, которые описывают насколько связаны два явления:

• Виды (столбцы) и экологические условия (сайты)
• Фенотипические классы (cтолбцы) и генотипы (строки)
• Фенотипы (столбцы) и географические выделы (строки)

Для примера используем задачку из классической генетики, где анализируется цвет и форма горошин

Пусть наши данные выглядят вот так

Представим эти данные в виде матрицы

peas <- matrix(c(99, 42, 29, 13), byrow = T, ncol = 2)
peas

##      [,1] [,2]
## [1,]   99   42
## [2,]   29   13

В ячейках содержится наблюдаемая частота - O (от “observed”) события X (цвет горошины) при условии события Y (форма горошины).

Мы хотим проверить есть ли связь между этими двумя явлениями.
То есть, на языке генетиков, сцеплены ли гены X и Y.

Мы хотим выяснить, являются ли цвет и форма семян независимыми признаками.

Если эти два признака независимы, то наша нулевая модель выглядит так \[ M_0: Ft = Ft\cdot p_{y,w} + Ft \cdot p_{y,s} + Ft \cdot p_{g,w} + Ft \cdot p_{g,s} \]

где
\(Ft\) - общая численность
\(p_{y,w}\) - Вероятность встретить желтых морщинистых
\(p_{y,s}\) - Вероятность встретить желтых гладких
\(p_{g,w}\) - Вероятность встретить зеленых морщинистых
\(p_{g,s}\) - Вероятность встретить желтых гладких

Для вычисления ожидаемых частот - E (от “expected”), при условии справедливости нулевой модели, нам нужно сделать следующее:

Найти общую численность семян \(Ft\) - это просто число, равное сумме всех ячеек.

Найти общие численности гладких и морщинистых семян \(f_i\) - это вектор из двух чисел, равных суммам строк (маргинальная сумма строк)

Найти общие численности желтых и зеленых семян \(f_j\) - это вектор из двух чисел, равных суммам столбцов (маргинальная сумма столбцов)

Ft <- sum(peas)

f_i <- apply(peas, 1, FUN = sum)

f_j <- apply(peas, 2, FUN = sum)

Оценка вероятности быть гладким: \(p_s = \frac{141} {183}\)

Оценка вероятности быть морщинистым: \(p_w = \frac{42} {183}\)

Оценка вероятности быть желтым: \(p_y = \frac{128} {183}\)

Оценка вероятности быть зеленым: \(p_g = \frac{55} {183}\)

В векторном виде

p_i <- f_i / Ft #Вектор вероятностей для формы
p_j <- f_j / Ft #Вектор вероятностей для цвета

Вероятности сочетаний признаков

Но! в этой таблице легко узнать результат матричного произведения двух векторов

\[ \textbf{a} \cdot \textbf{b} = \begin{pmatrix} p_y \\ p_g \end{pmatrix} \times \begin{pmatrix} p_s & p_w \end{pmatrix} \]

Матрица вероятностей сочетаний

q <- p_i %*% t(p_j) 
q

##       [,1]  [,2]
## [1,] 0.539 0.232
## [2,] 0.161 0.069

Матрица ожидаемых частот

round(q * Ft, 1)

##      [,1] [,2]
## [1,] 98.6 42.4
## [2,] 29.4 12.6

Из курса статистики мы помним критерий \(\chi^2\)

\[ \chi_{total}^2 = \sum {\frac{(O_{ij} - E_{ij})^2} {E_{ij}}} \]

Вклад каждой ячейки в формирование общего \(\chi_{total}^2\)

\[ \chi_{ij} = \frac{ (O_{ij} - E_{ij})} {\sqrt E_{ij} } = \sqrt{Ft} \left[\frac {p_{ij} - p_ip_j}{\sqrt{p_ip_j}}\right] \]

\(Ft\) - Сумма всех частот

Здесь \(\frac{(O_{ij} - E_{ij})} {\sqrt{E_{ij}}}\) - стандартизованный остаток от нулевой модели.

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

E <- (p_i %*% t(p_j) * Ft)
O <- peas

sum((O-E)^2/E)

## [1] 0.0209

Сравним с результатами специализированной функции

chisq.test(x = O, p = q, correct = F)

## 
##  Pearson's Chi-squared test
## 
## data:  O
## X-squared = 0.02, df = 1, p-value = 0.9

Чем больше значение \(\chi^2\) тем сильнее данные отклоняются от нулевой модели, тем выше сопряженность между двумя явлениями. В случае с горохом общая сопряженность очень низкая (закон Менделя!)

Струкутра данных

Таблица данных содержит данные по численности 35 видов в 70 пробах.

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

\[ \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 отличается от матрицы, содержащей \(\chi_{ij}\), только тем, что каждое ее значение разделено на \(\sqrt{Ft}\). То есть это будут стандартизованные значения (что очень полезно, так как собственные значения будут меньше или равны 1).

В этой матрице \(p_i\) и \(p_j\) маргинальные вероятности, то есть

\(p_i\) - Вероятность встретить особь данного вида
\(p_j\) - Вероятность встретить любую особь в данной пробе

\(p_{ij}\) - вероятность встретить особь данного вида в данной пробе

Частоты и вероятности

f_ij <- mite #Частота встречи данного вида в данной пробе, то есть это первичные даные!

p_ij <- mite/Ft #вероятность встречи данного вида в данной пробе

q <- p_i %*% t(p_j) #вероятность встретить особь в данной пробе. 

Ft <- sum(mite) #Общее количество найденных животных

f_i <- apply(mite, MARGIN = 1, FUN = sum) #Общее количество особей в каждой пробе

p_i <- f_i/Ft #Вектор вероятностей встретить какую-либо особь в данной пробе

f_j <- apply(mite, MARGIN = 2, FUN = sum) #Общее количество особей в каждом виде

p_j <- f_j/Ft #Вектор вероятностей встретить особь данного вида

Матрица относительных частот p_ij для биолога имеет определенный смысл - это оценки вероятностей встречи особи данного вида в данной пробе.

НО! Терминология корреспондентного анализа берет свое начало в прикладной математике, во многом разработанной физиками. В терминологии корреспондентного анализа, суммы по строкам и столбцам в матрице относительных частот называются массой строки и массой столбца.

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

Момент инерции - интеграл элемента массы умноженной на квадрат расстояния до центра масс.

• Простим - это математикам…

Ожидаемые частоты для нулевой модели, то есть при условии, что все пробы и все виды независимы

E <- (p_i %*% t(p_j) * Ft) # Ожидаемые частоты

O <- mite # Наблюдаемые частоты (первичные данные)

Chi2 <- sum((O-E)^2/E)
Chi2

## [1] 16621

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

Вычислим величину, которую называют нормированным хи-квадратом, или “инерцией”

\[ Inertia = \frac{\chi^2}{Ft} \]

Она показывает удельный вклад каждой особи в формирование общего отклонения от нулевой модели

Inertia <- Chi2/Ft
Inertia

## [1] 1.7

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

#Матрица, вычисленная через вероятности
Q1 <- (p_ij - p_i %*% t(p_j))/sqrt(p_i %*% t(p_j)) 

#Та же матрица, вычисленная через частоты
Q <- (f_ij*Ft - f_i %*% t(f_j))/(Ft*sqrt(f_i %*% t(f_j))) 

Q <- as.matrix(Q)

sum(Q^2)

## [1] 1.7

Это та же самая \(Inertia\)!

Инерция - это мера отклонения от нулевой модели!

Чем больше инерция, тем сильнее взаимосвязаны изучаемые явления.

Матрица \(\textbf{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

dim(U) 

## [1] 70 35

dim(D)

## [1] 35 35

dim(V)

## [1] 35 35

Qsvd <- U %*% D %*% t(V) #матрица "восстановленная" из "вспомогательных" матриц 

round(sum(Q - Qsvd)) #разность между исходной и "восстановленной" матрицами

## [1] 0

Матрица исходная и матрица восстановленная совпадают.
Все работает!

Однако, разложив исходную матрицу на три “вспомогательных”, мы еще не достигли результата - нам необходимо произвести ординацию в ортогональных осях. Это можно сделать только с помощью собственных чисел и собственных векторов (eigen values decomposition).

Этот анализ можно произвести только на основе квадратной матрицы ковариации.

Как найти матрицу ковариации для матрицы вкладов в формирование общей сопряженности?

Вспомним: Собственные числа и собственные векторы мы искали для матрицы ассоциации (ковариации).

Эта матрица \(\textbf{A} = \textbf{X}' \textbf{X}\)

Где \(\textbf{X}\) - это центрированная матрица

Но! Центрованная матрица, по сути, и есть матрица отклонений от нулевой модели.

Нулевая модель в случае эвклдового простраства - Значения признака у всех объектов равны среднему.

В случае матрицы сопряжения \(\textbf{Q}\) - это тоже матрица отклонений от нулевой модели, поэтому

\[ \textbf {A} = \textbf{Q}'\textbf{Q} \]

Матрица A - квадратная симметричная матрица

A <- t(Q) %*% Q

eig_values <- eigen(A)$values #Собственные числа матрицы A
eig_vectors <- eigen(A)$vectors #Матрица собственных векторов для матрицы A 

qplot(eig_values, diag(D)) + ggtitle("Квадратичная зависимость!" )

Квадраты сингулярных чисел матрицы Q равны собственным значениям матрицы ассоциации

qplot(eig_values, diag(D)^2 ) + geom_abline()

Обратите внимание на две особенности

1. Собственные значения для матрицы ассоциации A = Q’Q равны квадратам сингулярных чисел для матрицы Q
   
2. Последнее собственное значение матрицы A равно нулю! Таково ее свойство (это связано с тем, что частоты в матрице Q центрированы). Поэтому в дальнейшем анализе будет использовано на одно собственное значение меньше.

Вспомним, что мы вычислили показатель \(Inertia\), как сумму квадратов матрицы Q, или, что тоже самое, \(Inertia\) - это удельный \(\chi^2\). Из этого пока не следует, что \(Inertia\) показывает общую изменчивость в системе.

Вспомним, что общая дисперсия в системе - это сумма собственных значений. Ранее мы вычислили вектор eig_values, содержащий собственные значения матрицы ассоциации A = Q’Q. Найдем сумму всех собственных значений, то есть суммарную дисперсию.

sum(eig_values)

## [1] 1.7

Это в точности значение \(Inertia\)!

• Главные оси независимы друг от друга (перпендикулярны)
• Каждая последующая объясняет меньше общей инерции (общего \(\chi^2\))
• Всего осей может быть не больше чем минимальное из этих значений: (число строк - 1), (число столбцов - 1). Число осей зависит от меньшего из измерений. Вспомним! Строки и столбцы в матрице, которую мы изучаем равноправны.
  
• Первая ось - переменные, которые объясняют максимум зависимости строк от столбцов (значения которых сильнее всего отличаются от ожидаемых для данных объектов)
• Результаты изображаются в виде точечных графиков, похожих на биплоты.

В PCA мера расстояния между точками n-мерном пространстве – эвклидово расстояние

\[ D = \sqrt{\sum(y_{i,j} - y_{i,k})^2} \]

В CA мера расстояния между точками в n-мерном пространстве – расстояние \(\chi^2\)

\[ \chi^2 = \sqrt{ \sum {\frac{1}{f_{+j}/Ft}} \left( \frac{f_{1j}}{f_{1+}} - \frac{f_{2j}}{f_{2+}} \right)^2} \]

При ординации в CA должны сохраняться расстояния \(\chi^2\) расстояния

Но строки и столбцы равноправны! Следовательно при совместной ординации надо выбрать способ центрирования: либо по столбцам (scaling = 2, или scaling = “species” ), либо по строкам (scaling = 1, или scaling = “sites”)

scaling 1 сохраняет \(\chi^2\) расстояния между объектами (строками)

scaling 2 сохраняет \(\chi^2\) расстояния между признаками (столбцами)

Становится виден отчетливый градиент!

Проведите корреспондентный анализ для данных по птицам Австралии.

bird_ca <- cca(birds[ , -c(1, 2)], scale = TRUE)
# summary(bird_pca)
screeplot(bird_ca,  bstick = TRUE) # график собственных чисел

# install.packages("remotes")
# remotes::install_github("gavinsimpson/ggvegan")

library(ggvegan)
autoplot(bird_ca, scaling = "sites") 

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

• 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 10.21.13).
• Quinn, G.G.P., Keough, M.J., 2002. Experimental design and data analysis for biologists. Cambridge University Press.
• Zuur, A.F., Ieno, E.N., Smith, G.M., 2007. Analysing ecological data. Springer.