Como @amoeba mencionó en los comentarios, PCA solo analizará un conjunto de datos y le mostrará los principales patrones (lineales) de variación en esas variables, las correlaciones o covarianzas entre esas variables y las relaciones entre las muestras (las filas ) en su conjunto de datos.

$\mathrm{Ca^{2+}}$

0,5 \times pag H + 1.4 \times C {una}^{2 +} + 0.1 \times T o t una l C una r si o norte

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

y el segundo componente

2.7 \times pag H + 0,3 \times C {una}^{2 +} - 5.6 \times T o t una l C una r si o norte

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

Estos componentes se pueden seleccionar libremente de las variables medidas, y las que se eligen son aquellas que explican secuencialmente la mayor cantidad de variación en el conjunto de datos, y que cada combinación lineal es ortogonal (no correlacionada con) las otras.

En una ordenación restringida, tenemos dos conjuntos de datos, pero no somos libres de seleccionar las combinaciones lineales del primer conjunto de datos (los datos químicos del suelo anteriores) que queramos. En cambio, tenemos que seleccionar combinaciones lineales de las variables en el segundo conjunto de datos que mejor expliquen la variación en el primero. Además, en el caso de PCA, el único conjunto de datos es la matriz de respuesta y no hay predictores (se podría pensar que la respuesta se predice a sí misma). En el caso restringido, tenemos un conjunto de datos de respuesta que deseamos explicar con un conjunto de variables explicativas.

Aunque no ha explicado qué variables son la respuesta, normalmente se desea explicar la variación en la abundancia o composición de esas especies (es decir, las respuestas) utilizando las variables explicativas ambientales.

La versión restringida de PCA es una cosa llamada Análisis de Redundancia (RDA) en círculos ecológicos. Esto supone un modelo de respuesta lineal subyacente para la especie, que no es apropiado o solo apropiado si tiene gradientes cortos a lo largo de los cuales la especie responde.

Una alternativa a PCA es una cosa llamada análisis de correspondencia (CA). Esto no tiene restricciones, pero tiene un modelo de respuesta unimodal subyacente, que es algo más realista en términos de cómo las especies responden a lo largo de gradientes más largos. Tenga en cuenta también que CA modela abundancias relativas o composición , PCA modela las abundancias en bruto.

Existe una versión restringida de CA, conocida como análisis de correspondencia restringida o canónica (CCA), que no debe confundirse con un modelo estadístico más formal conocido como análisis de correlación canónica.

Tanto en RDA como en CCA, el objetivo es modelar la variación en la abundancia o composición de especies como una serie de combinaciones lineales de las variables explicativas.

De la descripción parece que su esposa quiere explicar la variación en la composición (o abundancia) de especies de milpiés en términos de las otras variables medidas.

Algunas palabras de advertencia; RDA y CCA son solo regresiones multivariadas; CCA es solo una regresión multivariada ponderada. Todo lo que has aprendido sobre la regresión se aplica, y también hay un par de otras trampas:

A medida que aumenta el número de variables explicativas, las restricciones se vuelven cada vez menos y realmente no está extrayendo componentes / ejes que explican la composición de las especies de manera óptima
con CCA, a medida que aumenta el número de factores explicativos, corre el riesgo de inducir un artefacto de una curva en la configuración de puntos en la gráfica de CCA.
La teoría subyacente en RDA y CCA está menos desarrollada que los métodos estadísticos más formales. Solo podemos elegir razonablemente qué variables explicativas seguir usando la selección por pasos (que no es ideal por todas las razones por las que no nos gusta como método de selección en regresión) y tenemos que usar pruebas de permutación para hacerlo.

entonces mi consejo es el mismo que con la regresión; piense con anticipación cuáles son sus hipótesis e incluya variables que reflejen esas hipótesis. No solo arroje todas las variables explicativas a la mezcla.

Ejemplo

Ordenación sin restricciones

PCA

Mostraré un ejemplo comparando PCA, CA y CCA usando el paquete vegano para R que ayudo a mantener y que está diseñado para adaptarse a este tipo de métodos de ordenación:

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 no estandariza la inercia, a diferencia de Canoco, por lo que la varianza total es 1826 y los valores propios están en esas mismas unidades y suman 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

También vemos que el primer valor propio es aproximadamente la mitad de la varianza y con los dos primeros ejes hemos explicado ~ 80% de la varianza total

> 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

Se puede extraer un biplot de los puntajes de las muestras y especies en los dos primeros componentes principales.

> plot(pcfit)

ingrese la descripción de la imagen aquí

Hay dos problemas aquí

La ordenación está dominada esencialmente por tres especies, estas especies se encuentran más alejadas del origen, ya que estos son los taxones más abundantes en el conjunto de datos.
Hay un fuerte arco de curva en la ordenación, que sugiere un gradiente único largo o dominante que se ha desglosado en los dos componentes principales principales para mantener las propiedades métricas de la ordenación.

California

Una AC podría ayudar con estos dos puntos, ya que maneja mejor el gradiente largo debido al modelo de respuesta unimodal, y modela la composición relativa de especies, no de abundancias en bruto.

El código vegano / R para hacer esto es similar al código PCA utilizado anteriormente

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)

Aquí explicamos aproximadamente el 40% de la variación entre sitios en su composición relativa

> 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

La parcela conjunta de las especies y los puntajes del sitio ahora está menos dominada por algunas especies

> plot(cafit)

ingrese la descripción de la imagen aquí

El tipo de PCA o CA que elija debe estar determinado por las preguntas que desee hacer sobre los datos. Por lo general, con los datos de especies, a menudo estamos interesados en la diferencia en el conjunto de especies, por lo que CA es una opción popular. Si tenemos un conjunto de datos de variables ambientales, por ejemplo, química del agua o del suelo, no esperaríamos que respondan de manera unimodal a lo largo de los gradientes, por lo que CA sería inapropiado y PCA (de una matriz de correlación, usando scale = TRUEen la rda()llamada) sería más apropiado.

Ordenación restringida; CCA

Ahora, si tenemos un segundo conjunto de datos que deseamos usar para explicar patrones en el primer conjunto de datos de especies, debemos usar una ordenación restringida. A menudo, la elección aquí es CCA, pero RDA es una alternativa, al igual que RDA después de la transformación de los datos para permitirle manejar mejor los datos de especies.

data(varechem)                          # load explanatory example data

Reutilizamos la cca()función pero suministramos dos marcos de datos ( Xpara especies y Ypara variables explicativas / predictoras) o una fórmula modelo que enumera la forma del modelo que deseamos ajustar.

Para incluir todas las variables que podríamos usar varechem ~ ., data = varechemcomo fórmula para incluir todas las variables, pero como dije anteriormente, esta no es una buena idea en general

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

El triplot de la ordenación anterior se produce utilizando el plot()método

> plot(ccafit)

ingrese la descripción de la imagen aquí

Por supuesto, ahora la tarea es determinar cuál de esas variables es realmente importante. También tenga en cuenta que hemos explicado aproximadamente 2/3 de la varianza de las especies utilizando solo 13 variables. Uno de los problemas de usar todas las variables en esta ordenación es que hemos creado una configuración arqueada en los puntajes de muestra y especie, que es puramente un artefacto de usar demasiadas variables correlacionadas.

Si desea saber más sobre esto, consulte la documentación vegana o un buen libro sobre análisis de datos ecológicos multivariados.

Relación con regresión

Es más simple ilustrar el enlace con RDA, pero CCA es igual, excepto que todo involucra sumas marginales de tabla de dos vías de fila y columna como ponderaciones.

En esencia, la RDA es equivalente a la aplicación de PCA a una matriz de valores ajustados a partir de una regresión lineal múltiple ajustada a los valores de cada especie (respuesta) (abundancias, por ejemplo) con predictores dados por la matriz de variables explicativas.

En R podemos hacer esto como

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

Los valores propios de estos dos enfoques son iguales:

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

Por alguna razón, no puedo hacer que los puntajes del eje (cargas) coincidan, pero invariablemente estos se escalan (o no), así que necesito analizar exactamente cómo se están haciendo aquí.

No hacemos la RDA a través de rda()lo que mostré con lm()etc., pero usamos una descomposición QR para la parte del modelo lineal y luego SVD para la parte de PCA. Pero los pasos esenciales son los mismos.

Gavin Simpson
fuente

+ 1.3

$+1.3$

- 5.6

$-5.6$

m

$m$

n

$n$

n \times m

$n \times m$

@amoeba Perdón por el retraso, pero agregué una sección a mi respuesta para intentar mostrar el enlace con la regresión y cómo RDA puede verse como un PCA de valores ajustados de una serie de regresiones lineales, una por variable de respuesta.

Gavin Simpson

X β

$\mathbf{X}\beta$

β

$\beta$ fitted()

X β

$\mathbf{X}\beta$

El origen de la RDA se debe a Rao (1964), que es un documento estadístico, por lo que debería ser apropiado.

Gavin Simpson

¿Qué criterios usar para separar las variables en variables explicativas y respuestas para los métodos de ordenación en ecología?

Respuestas:

Ejemplo

Ordenación sin restricciones

PCA

California

Ordenación restringida; CCA

Relación con regresión