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,doza2 +, TotalCarbon, może się okazać, że był to pierwszy składnik
0,5 × p H + 1,4 × Cza2 ++ 0,1 × T o t a l C a r b o n
i drugi składnik
2,7 × p H + 0,3 × Cza2 +- 5,6 × T o t a l C a r b o n
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)
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)
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)
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.