Как упомянуто в комментариях @amoeba, PCA будет рассматривать только один набор данных и покажет основные (линейные) закономерности изменения этих переменных, корреляции или ковариации между этими переменными и отношения между выборками (строками). ) в вашем наборе данных.

То, что обычно делают с набором данных о видах и набором потенциальных объясняющих переменных, - это соответствие ограниченному порядку. В PCA главные компоненты, оси на биплоте PCA, выводятся как оптимальные линейные комбинации всех переменных. Если вы запускали это на наборе данных по химии почвы с переменным рН, , TotalCarbon, вы можете обнаружить , что первый компонент был $\mathrm{Ca^{2+}}$

0,5 \times п ЧАС + 1.4 \times С a^{2 +} + 0,1 \times T о T a L С a р б о N

$0.5 \times \mathrm{pH} + 1.4 \times \mathrm{Ca^{2+}} + 0.1 \times \mathrm{TotalCarbon}$

и второй компонент

2,7 \times п ЧАС + 0,3 \times С a^{2 +} - 5,6 \times T о T a L С a р б о N

$2.7 \times \mathrm{pH} + 0.3 \times \mathrm{Ca^{2+}} - 5.6 \times \mathrm{TotalCarbon}$

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

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

Хотя вы не объяснили, какие переменные являются ответом, обычно желательно объяснить изменение численности или состава этих видов (т. Е. Ответы) с помощью переменных, объясняющих окружающую среду.

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

Альтернативой PCA является то, что называется анализом соответствия (CA). Это не ограничено, но у него есть базовая модель унимодального ответа, которая несколько более реалистична с точки зрения того, как виды реагируют на более длинные градиенты. Также обратите внимание, что CA моделирует относительную численность или состав , PCA моделирует исходную численность.

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

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

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

Несколько слов предупреждения; RDA и CCA - это просто многовариантные регрессии; CCA - это просто взвешенная многомерная регрессия. Все, что вы узнали о регрессии, применимо, и есть несколько других ошибок:

когда вы увеличиваете количество объясняющей переменной, ограничения на самом деле становятся все меньше и меньше, и вы на самом деле не извлекаете компоненты / оси, которые оптимально объясняют видовой состав, и
с CCA, когда вы увеличиваете количество объясняющих факторов, вы рискуете вызвать артефакт кривой в конфигурации точек на графике CCA.
теория, лежащая в основе RDA и CCA, менее развита, чем более формальные статистические методы. Мы можем только разумно выбирать, какие пояснительные переменные продолжать использовать пошаговый отбор (что не является идеальным по всем причинам, по которым нам не нравится его как метод выбора в регрессии), и мы должны использовать тесты перестановки для этого.

поэтому мой совет такой же, как с регрессией; заранее продумайте, каковы ваши гипотезы, и включите переменные, которые отражают эти гипотезы. не просто бросайте все объясняющие переменные в микс.

пример

Неограниченное рукоположение

PCA

Я покажу пример, сравнивающий PCA, CA и CCA, используя пакет vegan для R, который я помогаю поддерживать и который разработан для соответствия этим методам ординации:

library("vegan")                        # load the package
data(varespec)                          # load example data

## PCA
pcfit <- rda(varespec)
## could add `scale = TRUE` if variables in different units
pcfit

> pcfit
Call: rda(X = varespec)

              Inertia Rank
Total            1826     
Unconstrained    1826   23
Inertia is variance 

Eigenvalues for unconstrained axes:
  PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
983.0 464.3 132.3  73.9  48.4  37.0  25.7  19.7 
(Showed only 8 of all 23 unconstrained eigenvalues)

В отличие от Canoco, веганский не стандартизирует инерцию, поэтому общая дисперсия равна 1826, а собственные значения находятся в тех же единицах и составляют 1826

> cumsum(eigenvals(pcfit))
      PC1       PC2       PC3       PC4       PC5       PC6       PC7       PC8 
 982.9788 1447.2829 1579.5334 1653.4670 1701.8853 1738.8947 1764.6209 1784.3265 
      PC9      PC10      PC11      PC12      PC13      PC14      PC15      PC16 
1796.6007 1807.0361 1816.3869 1819.1853 1821.5128 1822.9045 1824.1103 1824.9250 
     PC17      PC18      PC19      PC20      PC21      PC22      PC23 
1825.2563 1825.4429 1825.5495 1825.6131 1825.6383 1825.6548 1825.6594

Мы также видим, что первое собственное значение составляет примерно половину дисперсии, а с первыми двумя осями мы объяснили ~ 80% от общей дисперсии

> head(cumsum(eigenvals(pcfit)) / pcfit$tot.chi)
      PC1       PC2       PC3       PC4       PC5       PC6 
0.5384240 0.7927453 0.8651851 0.9056821 0.9322031 0.9524749

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

> plot(pcfit)

введите описание изображения здесь

Здесь есть два вопроса

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

Калифорния

CA может помочь с обеими этими точками, так как он лучше справляется с длинным градиентом благодаря модели унимодального отклика и моделирует относительный состав видов, а не сырые численности.

Код Vegan / R для этого похож на код PCA, использованный выше

cafit <- cca(varespec)
cafit

> cafit <- cca(varespec)
> cafit
Call: cca(X = varespec)

              Inertia Rank
Total           2.083     
Unconstrained   2.083   23
Inertia is mean squared contingency coefficient 

Eigenvalues for unconstrained axes:
   CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
0.5249 0.3568 0.2344 0.1955 0.1776 0.1216 0.1155 0.0889 
(Showed only 8 of all 23 unconstrained eigenvalues)

Здесь мы объясняем около 40% различий между сайтами в их относительном составе

> head(cumsum(eigenvals(cafit)) / cafit$tot.chi)
      CA1       CA2       CA3       CA4       CA5       CA6 
0.2519837 0.4232578 0.5357951 0.6296236 0.7148866 0.7732393

Совместный график оценок видов и участков в настоящее время менее доминирует среди нескольких видов.

> plot(cafit)

введите описание изображения здесь

Какой из PCA или CA вы выбираете, должно определяться вопросами, которые вы хотите задать относительно данных. Обычно с данными о видах нас больше интересует разница в наборе видов, поэтому СА является популярным выбором. Если у нас есть набор данных переменных среды, скажем, химия воды или почвы, мы не ожидаем, что они будут реагировать унимодальным образом вдоль градиентов, поэтому CA будет неуместным, а PCA (корреляционной матрицы, используемой scale = TRUEв rda()вызове) будет более подходящий.

Ограниченное рукоположение; CCA

Теперь, если у нас есть второй набор данных, который мы хотим использовать для объяснения закономерностей в наборе данных первого вида, мы должны использовать ограниченную ординацию. Часто здесь выбирается CCA, но RDA является альтернативой, как и RDA после преобразования данных, чтобы позволить им лучше обрабатывать данные видов.

data(varechem)                          # load explanatory example data

Мы повторно используем cca()функцию, но мы либо предоставляем два фрейма данных ( Xдля видов и Yдля объясняющих / предикторных переменных), либо формулу модели, перечисляющую форму модели, которую мы хотим подогнать.

Чтобы включить все переменные, которые мы могли бы использовать varechem ~ ., data = varechemв качестве формулы, чтобы включить все переменные - но, как я уже сказал выше, это не очень хорошая идея в целом

ccafit <- cca(varespec ~ ., data = varechem)

> ccafit
Call: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn +
Zn + Mo + Baresoil + Humdepth + pH, data = varechem)

              Inertia Proportion Rank
Total          2.0832     1.0000     
Constrained    1.4415     0.6920   14
Unconstrained  0.6417     0.3080    9
Inertia is mean squared contingency coefficient 

Eigenvalues for constrained axes:
  CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7   CCA8   CCA9  CCA10  CCA11 
0.4389 0.2918 0.1628 0.1421 0.1180 0.0890 0.0703 0.0584 0.0311 0.0133 0.0084 
 CCA12  CCA13  CCA14 
0.0065 0.0062 0.0047 

Eigenvalues for unconstrained axes:
    CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8     CA9 
0.19776 0.14193 0.10117 0.07079 0.05330 0.03330 0.01887 0.01510 0.00949

Триплот вышеуказанного рукоположения производится по plot()методу

> plot(ccafit)

введите описание изображения здесь

Конечно, теперь задача состоит в том, чтобы определить, какая из этих переменных действительно важна. Также обратите внимание, что мы объяснили около 2/3 дисперсии видов, используя только 13 переменных. Одна из проблем использования всех переменных в этом порядке состоит в том, что мы создали арочную конфигурацию в выборке и оценках видов, что является чисто артефактом использования слишком большого количества коррелированных переменных.

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

Связь с регрессией

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

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

В R мы можем сделать это как

## centre the responses
spp <- scale(data.matrix(varespec), center = TRUE, scale = FALSE)
## ...and the predictors
env <- as.data.frame(scale(varechem, center = TRUE, scale = FALSE))

## fit a linear model to each column (species) in spp.
## Suppress intercept as we've centred everything
fit <- lm(spp ~ . - 1, data = env)

## Collect fitted values for each species and do a PCA of that
## matrix
pclmfit <- prcomp(fitted(fit))

Собственные значения для этих двух подходов равны:

> (eig1 <- unclass(unname(eigenvals(pclmfit)[1:14])))
 [1] 820.1042107 399.2847431 102.5616781  47.6316940  26.8382218  24.0480875
 [7]  19.0643756  10.1669954   4.4287860   2.2720357   1.5353257   0.9255277
[13]   0.7155102   0.3118612
> (eig2 <- unclass(unname(eigenvals(rdafit, constrained = TRUE))))
 [1] 820.1042107 399.2847431 102.5616781  47.6316940  26.8382218  24.0480875
 [7]  19.0643756  10.1669954   4.4287860   2.2720357   1.5353257   0.9255277
[13]   0.7155102   0.3118612
> all.equal(eig1, eig2)
[1] TRUE

По какой-то причине я не могу добиться совпадения баллов по осям (нагрузок), но они неизменно масштабируются (или нет), поэтому мне нужно посмотреть, как именно это делается здесь.

Мы не выполняем RDA через, rda()как я показал с помощью lm()etc, но мы используем QR-разложение для части линейной модели, а затем SVD для части PCA. Но основные шаги одинаковы.

Гэвин Симпсон
источник

+ 1.3

$+1.3$

- 5.6

$-5.6$

Спасибо за предложение и продолжение - мне не пришло в голову сделать линейные комбинации примеров ортогональными, но я обновил их. По поводу 2) я высказал предположение, но учитывая, что есть ц. 12 000 видов многоножек, я подозреваю, что ответ здесь - наблюдения об изобилии

m

$m$ виды на каждом из

n

$n$ места отбора проб. В этом смысле RDA или CCA будут моделировать многомерную матрицу ответов измерения

n \times m

$n \times m$ , Я постараюсь разобраться с 4 позже, после того как я уложу детей спать.

Гэвин Симпсон

@amoeba Извините за задержку, но я добавил раздел к своему ответу, чтобы попытаться показать связь с регрессией и то, как RDA можно рассматривать как PCA подходящих значений из серии линейных регрессий, по одной на переменную ответа.

Гэвин Симпсон

@amoeba Мы делаем SVD из

X β

$\mathbf{X}\beta$ (установленные значения) не коэффициентов

β

$\beta$ по крайней мере вот что fitted()дает:

X β

$\mathbf{X}\beta$ , Следовательно, RDA часто называют регрессией пониженного ранга.

Гэвин Симпсон

Происхождение RDA связано с Рао (1964), который является статистическим документом, поэтому должен быть уместным.

Гэвин Симпсон

Какие критерии использовать для разделения переменных на объясняющие переменные и ответы для методов ординации в экологии?

Ответы:

пример

Неограниченное рукоположение

PCA

Калифорния

Ограниченное рукоположение; CCA

Связь с регрессией