J'ai différentes variables qui interagissent au sein d'une population. Fondamentalement, j'ai fait un inventaire des mille-pattes et mesuré d'autres valeurs du terrain, comme:

• L'espèce et la quantité de spécimens collectés
• Les différents environnements où se trouvent les animaux
• le pH
• Le pourcentage de matière organique
• la quantité de P, K, Mg, Ca, Mn, Fe, Zn, Cu
• Relation Ca + Mg / K

Fondamentalement, je voudrais utiliser l'ACP pour déterminer quelles variables déterminent la variabilité des échantillons et rendent la forêt (les environnements) différente; quelles variables dois-je utiliser pour les "variables" et lesquelles pour les "individus"?

Comme @amoeba l'a mentionné dans les commentaires, l'ACP ne regardera qu'un seul ensemble de données et vous montrera les principaux modèles (linéaires) de variation de ces variables, les corrélations ou covariances entre ces variables et les relations entre les échantillons (les lignes ) dans votre ensemble de données.

Ce que l'on fait normalement avec un ensemble de données sur les espèces et une suite de variables explicatives potentielles, c'est d'adapter une ordination contrainte. En PCA, les composantes principales, les axes sur le biplot PCA, sont dérivées comme des combinaisons linéaires optimales de toutes les variables. Si vous avez exécuté ceci sur un ensemble de données de chimie du sol avec des variables pH, , TotalCarbon, vous pourriez constater que le premier composant était$\mathrm{Ca^{2+}}$

et le deuxième composant

Ces composants sont librement sélectionnables parmi les variables mesurées, et celles qui sont choisies sont celles qui expliquent séquentiellement la plus grande variation dans l'ensemble de données, et que chaque combinaison linéaire est orthogonale (non corrélée avec) les autres.

Dans une ordination contrainte, nous avons deux ensembles de données, mais nous ne sommes pas libres de sélectionner les combinaisons linéaires du premier ensemble de données (les données de chimie du sol ci-dessus) que nous voulons. Au lieu de cela, nous devons sélectionner des combinaisons linéaires des variables dans le deuxième ensemble de données qui expliquent le mieux la variation dans le premier. De plus, dans le cas de l'ACP, le seul ensemble de données est la matrice de réponse et il n'y a pas de prédicteurs (vous pourriez penser que la réponse se prédisait elle-même). Dans le cas contraint, nous avons un ensemble de données de réponse que nous souhaitons expliquer avec un ensemble de variables explicatives.

Bien que vous n'ayez pas expliqué quelles variables sont la réponse, on souhaite normalement expliquer la variation de l'abondance ou de la composition de ces espèces (c.-à-d. Les réponses) en utilisant les variables explicatives environnementales.

La version contrainte de PCA est une chose appelée analyse de redondance (RDA) dans les cercles écologiques. Cela suppose un modèle de réponse linéaire sous-jacent pour l'espèce, qui n'est pas approprié ou uniquement approprié si vous avez de courts gradients le long desquels l'espèce répond.

Une alternative à l'ACP est une chose appelée analyse des correspondances (AC). Ce n'est pas contraint, mais il a un modèle de réponse unimodal sous-jacent, qui est un peu plus réaliste en termes de réponse des espèces le long de gradients plus longs. Notez également que l'AC modélise les abondances ou la composition relatives , l'ACP modélise les abondances brutes.

Il existe une version contrainte de l'AC, connue sous le nom d' analyse de correspondance contrainte ou canonique (CCA) - à ne pas confondre avec un modèle statistique plus formel connu sous le nom d'analyse de corrélation canonique.

Tant dans la RDA que dans la CCA, l'objectif est de modéliser la variation de l'abondance ou de la composition des espèces comme une série de combinaisons linéaires des variables explicatives.

D'après la description, il semble que votre femme veuille expliquer la variation de la composition (ou de l'abondance) des espèces de mille-pattes en fonction des autres variables mesurées.

Quelques mots d'avertissement; RDA et CCA ne sont que des régressions multivariées; L'ACC n'est qu'une régression multivariée pondérée. Tout ce que vous avez appris sur la régression s'applique, et il existe également quelques autres pièges:

• à mesure que vous augmentez le nombre de variables explicatives, les contraintes deviennent de moins en moins et vous n'extrayez pas vraiment les composants / axes qui expliquent la composition des espèces de manière optimale, et
• avec CCA, en augmentant le nombre de facteurs explicatifs, vous risquez d'induire un artefact d'une courbe dans la configuration des points du tracé CCA.
• la théorie qui sous-tend la RDA et la CCA est moins bien développée que les méthodes statistiques plus formelles. Nous ne pouvons que raisonnablement choisir quelles variables explicatives continuer à utiliser la sélection par étapes (ce qui n'est pas idéal pour toutes les raisons pour lesquelles nous ne l'aimons pas comme méthode de sélection en régression) et nous devons utiliser des tests de permutation pour le faire.

donc mon conseil est le même que pour la régression; réfléchissez à l'avance à vos hypothèses et incluez des variables qui reflètent ces hypothèses. Ne vous contentez pas de jeter toutes les variables explicatives dans le mélange.

Exemple

Ordination sans contrainte

PCA

Je vais montrer un exemple comparant PCA, CA et CCA en utilisant le package vegan pour R que j'aide à maintenir et qui est conçu pour s'adapter à ces types de méthodes d'ordination:

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 ne standardise pas l'inertie, contrairement à Canoco, donc la variance totale est de 1826 et les valeurs propres sont dans ces mêmes unités et totalisent 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

Nous voyons également que la première valeur propre est environ la moitié de la variance et avec les deux premiers axes, nous avons expliqué ~ 80% de la variance totale

> 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

Un biplot peut être tiré des scores des échantillons et des espèces sur les deux premières composantes principales

> plot(pcfit)

Il y a deux problèmes ici

1. L'ordination est essentiellement dominée par trois espèces - ces espèces sont les plus éloignées de l'origine - car ce sont les taxons les plus abondants dans l'ensemble de données
2. Il y a une forte arche de courbe dans l'ordination, suggérant un gradient unique long ou dominant qui a été décomposé en deux principales composantes principales pour maintenir les propriétés métriques de l'ordination.

Californie

Une AC pourrait aider avec ces deux points car elle gère mieux le long gradient en raison du modèle de réponse unimodale, et elle modélise la composition relative des espèces et non les abondances brutes.

Le code vegan / R pour ce faire est similaire au code PCA utilisé ci-dessus

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) 

Nous expliquons ici environ 40% de la variation entre les sites dans leur composition relative

> 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

L'intrigue conjointe des scores des espèces et des sites est désormais moins dominée par quelques espèces

> plot(cafit)

Le PCA ou l'AC que vous choisissez doit être déterminé par les questions que vous souhaitez poser aux données. Habituellement, avec les données sur les espèces, nous sommes plus souvent intéressés par les différences dans la suite d'espèces, donc l'AC est un choix populaire. Si nous avons un ensemble de données de variables environnementales, par exemple la chimie de l'eau ou du sol, nous ne nous attendrions pas à ce que celles-ci répondent de manière unimodale le long des gradients, donc l'AC serait inapproprié et l'ACP (d'une matrice de corrélation, utilisée scale = TRUEdans l' rda()appel) serait plus approprié.

Ordination contrainte; CCA

Maintenant, si nous avons un deuxième ensemble de données que nous souhaitons utiliser pour expliquer les modèles dans le premier ensemble de données sur les espèces, nous devons utiliser une ordination contrainte. Souvent, le choix ici est CCA, mais RDA est une alternative, tout comme RDA après transformation des données pour lui permettre de mieux gérer les données sur les espèces.

data(varechem)                          # load explanatory example data

Nous réutilisons la cca()fonction mais nous fournissons soit deux bases de données ( Xpour les espèces et Ypour les variables explicatives / prédictives), soit une formule de modèle répertoriant la forme du modèle que nous souhaitons adapter.

Pour inclure toutes les variables, nous pourrions utiliser varechem ~ ., data = varechemcomme formule pour inclure toutes les variables - mais comme je l'ai dit ci-dessus, ce n'est pas une bonne idée en général

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

Le triplot de l'ordination ci-dessus est produit en utilisant la plot()méthode

> plot(ccafit)

Bien sûr, la tâche consiste maintenant à déterminer laquelle de ces variables est réellement importante. Notez également que nous avons expliqué environ 2/3 de la variance des espèces en utilisant seulement 13 variables. l'un des problèmes de l'utilisation de toutes les variables dans cette ordination est que nous avons créé une configuration arquée dans les scores des échantillons et des espèces, ce qui est purement un artefact de l'utilisation d'un trop grand nombre de variables corrélées.

Si vous voulez en savoir plus à ce sujet, consultez la documentation végétalienne ou un bon livre sur l'analyse des données écologiques multivariées.

Relation avec la régression

Il est plus simple d'illustrer le lien avec RDA, mais l'ACC est identique, sauf que tout implique des sommes marginales de table bidirectionnelle de ligne et de colonne comme poids.

En son cœur, la RDA équivaut à l'application de l'ACP à une matrice de valeurs ajustées à partir d'une régression linéaire multiple ajustée à chaque espèce (réponse) valeurs (abondances, disons) avec des prédicteurs donnés par la matrice de variables explicatives.

Dans R, nous pouvons le faire comme

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

Les valeurs propres de ces deux approches sont égales:

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

Pour une raison quelconque, je ne parviens pas à faire correspondre les scores d'axe (chargements), mais invariablement, ils sont mis à l'échelle (ou non), je dois donc examiner exactement comment ceux-ci sont effectués ici.

Nous ne faisons pas le RDA via rda()comme je l'ai montré avec lm()etc, mais nous utilisons une décomposition QR pour la partie modèle linéaire puis SVD pour la partie PCA. Mais les étapes essentielles sont les mêmes.