Wie in den Kommentaren von @amoeba erwähnt, betrachtet PCA nur einen Datensatz und zeigt Ihnen die wichtigsten (linearen) Variationsmuster dieser Variablen, die Korrelationen oder Kovarianzen zwischen diesen Variablen und die Beziehungen zwischen Stichproben (den Zeilen) ) in Ihrem Datensatz.

Was man normalerweise mit einem Artendatensatz und einer Reihe möglicher erklärender Variablen macht, ist, sich einer eingeschränkten Ordination anzupassen. In PCA werden die Hauptkomponenten, die Achsen im PCA-Biplot, als optimale lineare Kombinationen aller Variablen abgeleitet. Wenn Sie dies mit einem Datensatz der Bodenchemie mit den Variablen pH, , TotalCarbon, durchgeführt haben, stellen Sie möglicherweise fest, dass die erste Komponente war $\mathrm{Ca^{2+}}$

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}$

und die zweite Komponente

2.7 \times p H + 0.3 \times C a^{2 +} - 5.6 \times T o t a l C a r b o n

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

Diese Komponenten können aus den gemessenen Variablen frei ausgewählt werden. Es werden diejenigen ausgewählt, die nacheinander die größte Variation im Datensatz erklären und bei denen jede lineare Kombination orthogonal (nicht korreliert mit den anderen) ist.

In einer eingeschränkten Ordination haben wir zwei Datensätze, aber wir können nicht frei wählen, welche linearen Kombinationen des ersten Datensatzes (die obigen chemischen Bodendaten) wir wollen. Stattdessen müssen wir lineare Kombinationen der Variablen im zweiten Datensatz auswählen, die die Variation im ersten am besten erklären. Im Fall von PCA ist der eine Datensatz die Antwortmatrix, und es gibt keine Prädiktoren (Sie können sich die Antwort als Vorhersage selbst vorstellen). Im eingeschränkten Fall haben wir einen Antwortdatensatz, den wir mit einem Satz erklärender Variablen erklären möchten.

Obwohl Sie nicht erklärt haben, welche Variablen die Antwort sind, möchten Sie normalerweise die Variation der Häufigkeit oder Zusammensetzung dieser Arten (dh der Antworten) anhand der erklärenden Umgebungsvariablen erklären.

Die eingeschränkte Version von PCA wird in ökologischen Kreisen als Redundanzanalyse (RDA) bezeichnet. Dies setzt ein zugrunde liegendes lineares Antwortmodell für die Art voraus, das entweder nicht oder nur dann geeignet ist, wenn Sie kurze Gradienten haben, entlang derer die Art reagiert.

Eine Alternative zu PCA ist die sogenannte Korrespondenzanalyse (CA). Dies ist nicht eingeschränkt, aber es liegt ein unimodales Antwortmodell zugrunde, das etwas realistischer ist, was die Reaktion von Arten entlang längerer Gradienten betrifft. Beachten Sie auch, dass CA relative Häufigkeiten oder Zusammensetzung modelliert, PCA modelliert die rohen Häufigkeiten.

Es gibt eine eingeschränkte Version von CA, die als eingeschränkte oder kanonische Korrespondenzanalyse (CCA) bezeichnet wird - nicht zu verwechseln mit einem formaleren statistischen Modell, das als kanonische Korrelationsanalyse bezeichnet wird.

Sowohl in RDA als auch in CCA besteht das Ziel darin, die Variation der Artenhäufigkeit oder -zusammensetzung als eine Reihe linearer Kombinationen der erklärenden Variablen zu modellieren.

Aus der Beschreibung geht hervor, dass Ihre Frau die Variation der Zusammensetzung (oder Häufigkeit) der Tausendfüßer-Arten anhand der anderen gemessenen Variablen erklären möchte.

Einige warnende Worte; RDA und CCA sind nur multivariate Regressionen. CCA ist nur eine gewichtete multivariate Regression. Alles, was Sie über Regression gelernt haben, gilt, und es gibt noch ein paar andere Fallstricke:

Wenn Sie die Anzahl der erklärenden Variablen erhöhen, werden die Einschränkungen tatsächlich immer geringer und Sie extrahieren nicht wirklich Komponenten / Achsen, die die Artenzusammensetzung optimal erklären, und
Wenn Sie mit CCA die Anzahl der erklärenden Faktoren erhöhen, besteht die Gefahr, dass ein Artefakt einer Kurve in die Konfiguration der Punkte im CCA-Diagramm einbezogen wird.
Die Theorie, die RDA und CCA zugrunde liegt, ist weniger gut entwickelt als formalere statistische Methoden. Wir können nur vernünftigerweise auswählen, welche erklärenden Variablen für die schrittweise Auswahl verwendet werden sollen (was aus all den Gründen, die wir als Auswahlmethode bei der Regression nicht mögen, nicht ideal ist), und wir müssen dazu Permutationstests verwenden.

Mein Rat ist also der gleiche wie bei der Regression. Überlegen Sie sich im Voraus, was Ihre Hypothesen sind, und schließen Sie Variablen ein, die diese Hypothesen widerspiegeln. Werfen Sie nicht einfach alle erklärenden Variablen in die Mischung.

Beispiel

Uneingeschränkte Ordination

PCA

Ich zeige ein Beispiel für den Vergleich von PCA, CA und CCA mit dem veganen Paket für R, das ich pflege und das für diese Ordnungsmethoden geeignet ist:

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)

vegan standardisiert die Trägheit im Gegensatz zu Canoco nicht, daher beträgt die Gesamtvarianz 1826 und die Eigenwerte sind in denselben Einheiten und summieren sich zu 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

Wir sehen auch, dass der erste Eigenwert ungefähr die Hälfte der Varianz beträgt und mit den ersten beiden Achsen ~ 80% der Gesamtvarianz erklärt haben

> 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

Aus den Bewertungen der Proben und Arten der ersten beiden Hauptkomponenten kann ein Biplot gezogen werden

> plot(pcfit)

Geben Sie hier die Bildbeschreibung ein

Hier gibt es zwei Probleme

Die Ordination wird im Wesentlichen von drei Arten dominiert - diese Arten liegen am weitesten vom Ursprung entfernt -, da dies die am häufigsten vorkommenden Taxa im Datensatz sind
Die Ordination weist einen starken Kurvenbogen auf, der auf einen langen oder dominanten einzelnen Gradienten hinweist, der in die beiden Hauptkomponenten unterteilt wurde, um die metrischen Eigenschaften der Ordination beizubehalten.

CA.

Eine CA kann bei beiden Punkten hilfreich sein, da sie aufgrund des unimodalen Antwortmodells einen langen Gradienten besser handhabt und die relative Zusammensetzung von Arten modelliert, die nicht roh vorkommen.

Der vegane / R- Code hierfür ähnelt dem oben verwendeten PCA-Code

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)

Hier erklären wir ungefähr 40% der Variation zwischen Standorten in ihrer relativen Zusammensetzung

> 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

Die gemeinsame Darstellung der Arten- und Standortwerte wird jetzt weniger von einigen wenigen Arten dominiert

> plot(cafit)

Geben Sie hier die Bildbeschreibung ein

Welche PCA oder CA Sie wählen, sollte durch die Fragen bestimmt werden, die Sie an die Daten stellen möchten. Normalerweise sind wir bei Artendaten häufiger an Unterschieden in der Artenreihe interessiert, daher ist CA eine beliebte Wahl. Wenn wir einen Datensatz von Umgebungsvariablen haben, z. B. Wasser- oder Bodenchemie, würden wir nicht erwarten, dass diese entlang von Gradienten unimodal reagieren, sodass CA unangemessen wäre und PCA (einer Korrelationsmatrix, die scale = TRUEim rda()Aufruf verwendet wird) angemessener.

Eingeschränkte Ordination; CCA

Wenn wir nun einen zweiten Datensatz haben, mit dem wir Muster im ersten Speziesdatensatz erklären möchten, müssen wir eine eingeschränkte Ordination verwenden. Oft ist hier CCA die Wahl, aber RDA ist eine Alternative, ebenso wie RDA nach der Transformation der Daten, damit sie besser mit Artendaten umgehen können.

data(varechem)                          # load explanatory example data

Wir verwenden die cca()Funktion wieder, geben jedoch entweder zwei Datenrahmen ( Xfür Arten und Yfür erklärende / prädiktive Variablen) oder eine Modellformel an, in der die Form des Modells aufgeführt ist, das angepasst werden soll.

Um alle Variablen einzuschließen, könnten wir varechem ~ ., data = varechemals Formel verwenden, um alle Variablen einzuschließen - aber wie ich oben sagte, ist dies im Allgemeinen keine gute Idee

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

Das Triplot der obigen Ordination wird unter Verwendung des plot()Verfahrens erzeugt

> plot(ccafit)

Geben Sie hier die Bildbeschreibung ein

Jetzt besteht die Aufgabe natürlich darin, herauszufinden, welche dieser Variablen tatsächlich wichtig ist. Beachten Sie auch, dass wir ungefähr 2/3 der Artenvarianz mit nur 13 Variablen erklärt haben. Eines der Probleme bei der Verwendung aller Variablen in dieser Ordination besteht darin, dass wir eine gewölbte Konfiguration in Stichproben- und Artenbewertungen erstellt haben, die lediglich ein Artefakt der Verwendung zu vieler korrelierter Variablen ist.

Wenn Sie mehr darüber erfahren möchten, lesen Sie die vegane Dokumentation oder ein gutes Buch über multivariate ökologische Datenanalysen.

Beziehung zur Regression

Es ist am einfachsten, den Zusammenhang mit RDA zu veranschaulichen, aber CCA ist genauso, außer dass alles Rand- und Spalten-Zwei-Wege-Tabellensummen als Gewichte umfasst.

Im Kern entspricht RDA der Anwendung von PCA auf eine Matrix angepasster Werte aus einer multiplen linearen Regression, die an die Werte jeder Art (Antwort) (z. B. Häufigkeiten) angepasst ist, wobei Prädiktoren durch die Matrix der erklärenden Variablen angegeben werden.

In R können wir dies als tun

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

Die Eigenwerte für diese beiden Ansätze sind gleich:

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

Aus irgendeinem Grund kann ich nicht erreichen, dass die Achsenwerte (Ladungen) übereinstimmen, aber diese sind immer skaliert (oder nicht), daher muss ich genau untersuchen, wie diese hier durchgeführt werden.

Wir machen die RDA nicht über, rda()wie ich mit lm()etc gezeigt habe, aber wir verwenden eine QR-Zerlegung für den linearen Modellteil und dann SVD für den PCA-Teil. Die wesentlichen Schritte sind jedoch dieselben.

Gavin Simpson
quelle

+1 und auf der Suche nach Fortsetzung! Mehrere Kommentare: (1) In Ihrem Beispiel ist PC1 nicht orthogonal zu PC2. Sie könnten zB auf ändern, um das zu beheben. (2) Es wäre wahrscheinlich sinnvoll, den Titel des OP so zu bearbeiten, dass er den Inhalt Ihrer Antwort widerspiegelt. Die aktuelle Ausgabe des Titels gehört mir, aber ich hatte keine Ahnung, wovon OP sprach. (3) Ist die "Antwort" normalerweise univariat oder multivariat? Klingt wie letzteres, aber wie ist zB Tausendfüßlerhäufigkeit multivariat? Überfluss an mehreren Arten? (4) Wie unterscheiden sich diese Methoden von der Regression? Können Sie einige mathematische Hinweise einfügen?

+ 1.3

$+1.3$

- 5.6

$-5.6$

Amöbe

Vielen Dank für den Vorschlag und das Follow-up - es ist mir nicht in den Sinn gekommen, die linearen Kombinationsbeispiele orthogonal zu machen, aber ich habe sie aktualisiert. Zu 2) habe ich eine Vermutung gemacht, aber da es c gibt. 12.000 Tausendfüßerarten, ich vermute, die Antwort hier sind Beobachtungen der Häufigkeit von Arten an jedem der Probenahmestellen. In diesem Sinne würde die RDA oder CCA eine multivariate Antwortmatrix mit der Dimension modellieren . Ich werde später versuchen, mit 4 fertig zu werden, nachdem ich die Kinder ins Bett gebracht habe.

m

$m$

n

$n$

n \times m

$n \times m$

Gavin Simpson

@amoeba Entschuldigen Sie die Verzögerung, aber ich habe meiner Antwort einen Abschnitt hinzugefügt, um zu versuchen, den Zusammenhang mit der Regression zu zeigen und wie RDA als PCA angepasster Werte aus einer Reihe linearer Regressionen angesehen werden kann, eine pro Antwortvariable.

Gavin Simpson

@amoeba Wir machen SVD von (angepasste Werte), nicht von den Koeffizienten , zumindest ergibt das: . Daher wird RDA oft als Regression mit reduziertem Rang bezeichnet.

X β

$\mathbf{X}\beta$

β

$\beta$ fitted()

X β

$\mathbf{X}\beta$

Gavin Simpson

Der Ursprung der RDA liegt bei Rao (1964) , einem statistischen Papier, das angemessen sein sollte.

Gavin Simpson

Welche Kriterien sind für die Trennung von Variablen in erklärende Variablen und Antworten für Ordnungsmethoden in der Ökologie zu verwenden?

Antworten:

Beispiel

Uneingeschränkte Ordination

PCA

CA.

Eingeschränkte Ordination; CCA

Beziehung zur Regression