Jak wspomniano w komentarzach @amoeba, PCA przyjrzy się tylko jednemu zestawowi danych i pokaże główne (liniowe) wzorce zmienności tych zmiennych, korelacje lub kowariancje między tymi zmiennymi oraz relacje między próbkami (wiersze ) w zestawie danych.

To, co zwykle robi się z zestawem danych o gatunkach i pakietem potencjalnych zmiennych objaśniających, to dopasowanie ograniczonego wyświęcenia. W PCA główne komponenty, osie na biplocie PCA, są wyprowadzane jako optymalne kombinacje liniowe wszystkich zmiennych. Jeśli użyłeś tego w zestawie danych chemii gleby ze zmiennymi pH, $\mathrm{Ca^{2+}}$ , TotalCarbon, może się okazać, że był to pierwszy składnik

0.5 \times p H + 1.4 \times C a^{2 +} + 0.1 \times T o t a l C a r b o n

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

i drugi składnik

2.7 \times p H. + 0,3 \times do {za}^{2) +} - 5.6 \times T. o t za l do za r b o n

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

Składniki te można dowolnie wybierać spośród mierzonych zmiennych, a wybierane są te, które wyjaśniają sekwencyjnie największą ilość zmian w zbiorze danych oraz że każda kombinacja liniowa jest ortogonalna (nieskorelowana) z innymi.

W ograniczonym święceniu mamy dwa zestawy danych, ale nie jesteśmy w stanie wybrać dowolnych liniowych kombinacji pierwszego zestawu danych (dane dotyczące chemii gleby powyżej), jakie chcemy. Zamiast tego musimy wybrać liniowe kombinacje zmiennych w drugim zestawie danych, które najlepiej wyjaśniają zmienność w pierwszym. Ponadto w przypadku PCA jednym zestawem danych jest matryca odpowiedzi i nie ma predyktorów (można by pomyśleć o odpowiedzi jako o przewidywaniu). W przypadku ograniczonym mamy zestaw danych odpowiedzi, który chcemy wyjaśnić za pomocą zestawu zmiennych objaśniających.

Chociaż nie wyjaśniłeś, które zmienne są odpowiedzią, zwykle chce się wyjaśnić zmienność liczebności lub składu tych gatunków (tj. Odpowiedzi) za pomocą zmiennych objaśniających środowisko.

Ograniczona wersja PCA jest nazywana analizą redundancji (RDA) w kręgach ekologicznych. Zakłada to podstawowy liniowy model reakcji dla gatunku, który jest albo nieodpowiedni, albo odpowiedni tylko, jeśli masz krótkie gradienty, wzdłuż których gatunek reaguje.

Alternatywą dla PCA jest tzw. Analiza korespondencji (CA). Jest to nieograniczone, ale ma bazowy model reakcji jednomodalnej, który jest nieco bardziej realistyczny pod względem reakcji gatunków na dłuższe gradienty. Należy również zauważyć, że CA modeluje względne obfitości lub skład , PCA modeluje surowe obfitości.

Istnieje ograniczona wersja CA, znana jako ograniczona lub kanoniczna analiza korespondencji (CCA) - której nie należy mylić z bardziej formalnym modelem statystycznym zwanym analizą korelacji kanonicznej.

Zarówno w RDA, jak i CCA celem jest modelowanie zmienności liczebności i składu gatunkowego jako serii liniowych kombinacji zmiennych objaśniających.

Z opisu brzmi to tak, jakby twoja żona chciała wyjaśnić zmienność składu gatunków krocionogów (lub liczebność) pod względem innych mierzonych zmiennych.

Kilka słów ostrzeżenia; RDA i CCA to tylko regresje wielowymiarowe; CCA to tylko ważona regresja wielowymiarowa. Wszystko, czego nauczyłeś się o regresji, ma zastosowanie, a także kilka innych błędów:

wraz ze wzrostem liczby zmiennych objaśniających ograniczenia stają się coraz mniejsze i tak naprawdę nie wyodrębniasz komponentów / osi, które optymalnie wyjaśniają skład gatunkowy, i
w przypadku CCA, gdy zwiększasz liczbę czynników objaśniających, ryzykujesz wywołaniem artefaktu krzywej do konfiguracji punktów na wykresie CCA.
teoria leżąca u podstaw RDA i CCA jest mniej rozwinięta niż bardziej formalne metody statystyczne. Możemy tylko rozsądnie wybierać zmienne objaśniające, aby nadal używać selekcji stopniowej (co nie jest idealne z wszystkich powodów, dla których nie lubimy jej jako metody selekcji w regresji) i musimy do tego celu użyć testów permutacyjnych.

więc moja rada jest taka sama jak w przypadku regresji; pomyśl z góry, jakie są twoje hipotezy i uwzględnij zmienne, które odzwierciedlają te hipotezy. Nie wrzucaj do mieszanki wszystkich zmiennych objaśniających.

Przykład

Święcenia nieograniczone

PCA

Pokażę przykład porównujący PCA, CA i CCA przy użyciu pakietu wegańskiego dla R, który pomagam utrzymać i który jest zaprojektowany, aby pasował do tego rodzaju metod święceń:

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)

wegańskie nie standaryzuje bezwładności, w przeciwieństwie do Canoco, więc całkowita wariancja wynosi 1826, a wartości własne są w tych samych jednostkach i sumują się do 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

Widzimy również, że pierwsza wartość własna stanowi około połowy wariancji, a przy dwóch pierwszych osiach wyjaśniliśmy ~ 80% całkowitej wariancji

> 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

Dwupłat można wyciągnąć z wyników próbek i gatunków na pierwszych dwóch głównych składnikach

> plot(pcfit)

wprowadź opis zdjęcia tutaj

Istnieją tutaj dwa problemy

Wyświęcenie jest zasadniczo zdominowane przez trzy gatunki - gatunki te znajdują się najdalej od pochodzenia - ponieważ są to najliczniejsze taksony w zbiorze danych
W wyświęceniu występuje silny łuk krzywej, sugerujący długi lub dominujący pojedynczy gradient, który został podzielony na dwa główne główne składniki w celu zachowania metrycznych właściwości wyświęcenia.

CA

CA może pomóc w obu tych punktach, ponieważ lepiej radzi sobie z długim gradientem ze względu na model reakcji jednomodalnej i modeluje względny skład gatunków, a nie surowe liczebności.

Kod wegański / R do tego celu jest podobny do kodu PCA użytego powyżej

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)

Tutaj wyjaśniamy około 40% zmienności między stronami w ich względnym składzie

> 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

Wspólna fabuła gatunków i oceny terenów jest obecnie mniej zdominowana przez kilka gatunków

> plot(cafit)

wprowadź opis zdjęcia tutaj

Wybór PCA lub urzędu certyfikacji powinien zależeć od pytań, które chcesz zadać na temat danych. Zwykle w przypadku danych gatunków częściej interesuje nas różnica w zestawie gatunków, więc CA jest popularnym wyborem. Gdybyśmy mieli zestaw danych zmiennych środowiskowych, powiedzmy chemii wody lub gleby, nie spodziewalibyśmy się, że reagują one w sposób jednomodalny wzdłuż gradientów, więc CA byłby niewłaściwy, a PCA (macierzy korelacji używanej scale = TRUEw rda()wywołaniu) byłby bardziej odpowiedni.

Ograniczone święcenia; CCA

Teraz, jeśli mamy drugi zestaw danych, który chcemy wykorzystać do wyjaśnienia wzorców w pierwszym zbiorze danych o gatunkach, musimy zastosować ograniczenie święceń. Często wybierany jest tutaj CCA, ale RDA jest alternatywą, podobnie jak RDA po transformacji danych, aby umożliwić lepsze przetwarzanie danych dotyczących gatunków.

data(varechem)                          # load explanatory example data

Ponownie używamy tej cca()funkcji, ale albo dostarczamy dwie ramki danych ( Xdla gatunków i Ydla zmiennych objaśniających / predykcyjnych) lub formułę modelu zawierającą formę modelu, który chcemy dopasować.

Aby uwzględnić wszystkie zmienne, moglibyśmy użyć varechem ~ ., data = varechemformuły do uwzględnienia wszystkich zmiennych - ale jak powiedziałem powyżej, ogólnie nie jest to dobry pomysł

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

Tryplot powyższego wyświęcenia jest tworzony przy użyciu plot()metody

> plot(ccafit)

wprowadź opis zdjęcia tutaj

Oczywiście teraz zadaniem jest ustalenie, która z tych zmiennych jest rzeczywiście ważna. Zauważ też, że wyjaśniliśmy około 2/3 wariancji gatunkowych, używając tylko 13 zmiennych. jednym z problemów związanych ze stosowaniem wszystkich zmiennych w tym uporządkowaniu jest to, że stworzyliśmy łukową konfigurację w próbkach i ocenach gatunków, co jest czystym artefaktem użycia zbyt wielu skorelowanych zmiennych.

Jeśli chcesz dowiedzieć się więcej na ten temat, sprawdź wegańską dokumentację lub dobrą książkę na temat wielowymiarowej analizy danych ekologicznych.

Związek z regresją

Najłatwiej jest zilustrować związek z RDA, ale CCA jest dokładnie taki sam, z tym wyjątkiem, że wszystko obejmuje obustronne sumy krańcowe wierszy i kolumn jako wagi.

W istocie RDA jest równoważne zastosowaniu PCA do macierzy dopasowanych wartości z wielokrotnej regresji liniowej dopasowanej do wartości każdego gatunku (odpowiedzi) (powiedzmy obfitości) z predyktorami podanymi przez macierz zmiennych objaśniających.

W R możemy to zrobić jako

## 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))

Wartości własne dla tych dwóch podejść są równe:

> (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

Z jakiegoś powodu nie mogę dopasować wyników osi (obciążeń), ale niezmiennie są one skalowane (lub nie), więc muszę przyjrzeć się dokładnie, jak się to robi.

Nie wykonujemy RDA poprzez, rda()jak pokazałem z lm()itp., Ale używamy rozkładu QR dla części modelu liniowego, a następnie SVD dla części PCA. Ale niezbędne kroki są takie same.

Gavin Simpson
źródło

+1 i szukam kontynuacji! Kilka komentarzy: (1) w twoim przykładzie PC1 nie jest prostopadły do PC2; możesz np. zmienić

+ 1.3

$+1.3$ do

- 5.6

$-5.6$ naprawić to. (2) Prawdopodobnie miałoby sens edytowanie tytułu PO, aby odzwierciedlić treść Twojej odpowiedzi; obecna edycja tego tytułu jest moja, ale nie miałem pojęcia, o czym mówi OP. (3) Czy „odpowiedź” jest zwykle jedno- lub wielowymiarowa? Brzmi jak ta ostatnia, ale jak np. Wielowariantowa obfitość krocionogów? Obfitość kilku gatunków? (4) Czym różnią się te metody od regresji? Czy możesz podać jakieś matematyczne wskazówki?

ameba

Dzięki za sugestię i kontynuację - nie przyszło mi do głowy, aby przykłady kombinacji liniowej były ortogonalne, ale je zaktualizowałem. Re 2), założyłem domniemanie, ale biorąc pod uwagę, że istnieją c. 12 000 gatunków krocionogów, podejrzewam, że odpowiedzią są obserwacje obfitości

m

$m$ gatunki na każdym z

n

$n$ miejsca pobierania próbek. W tym sensie RDA lub CCA modeluje wielowymiarową macierz odpowiedzi w wymiarze

n \times m

$n \times m$ . Spróbuję poradzić sobie z 4 później, kiedy położę dzieci do łóżka.

Gavin Simpson,

@amoeba Przepraszam za opóźnienie, ale dodałem sekcję do mojej odpowiedzi, aby spróbować pokazać link z regresją i jak RDA można postrzegać jako PCA dopasowanych wartości z serii regresji liniowych, po jednej na zmienną odpowiedzi.

Gavin Simpson,

@amoeba Robimy SVD z

X β

$\mathbf{X}\beta$ (dopasowane wartości) nie współczynników

β

$\beta$ , przynajmniej to fitted()daje:

X β

$\mathbf{X}\beta$ . Dlatego RDA jest często nazywane regresją o zmniejszonym stopniu rangi.

Gavin Simpson,

Pochodzenie RDA wynika z Rao (1964), który jest opracowaniem statystycznym, więc powinien być odpowiedni.

Gavin Simpson,

Jakie kryteria zastosować do podziału zmiennych na zmienne objaśniające i odpowiedzi na metody święceń w ekologii?

Odpowiedzi:

Przykład

Święcenia nieograniczone

PCA

CA

Ograniczone święcenia; CCA

Związek z regresją