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

Question

Tengo diferentes variables que interactúan dentro de una población. Básicamente he estado haciendo un inventario de milpiés y midiendo algunos otros valores del terreno, como:

• Las especies y la cantidad de especímenes recogidos
• Los diferentes entornos en los que se encuentran los animales
• el pH
• El porcentaje de materia orgánica
• la cantidad de P, K, Mg, Ca, Mn, Fe, Zn, Cu
• Relación Ca + Mg/K

Básicamente, me gustaría utilizar el ACP para determinar qué variables impulsan la variabilidad de las muestras y hacen que el bosque (entornos) sea diferente; ¿qué variables debo utilizar para las "variables" y cuáles para los "individuos"?

Answer 1

Como mencionó @amoeba en los comentarios, el PCA sólo mirará 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.

Lo que se hace normalmente con un conjunto de datos de especies y un conjunto de posibles variables explicativas es ajustar una ordenación restringida. En PCA, los componentes principales, los ejes en el biplot de PCA, se derivan como combinaciones lineales óptimas de todas las variables. Si se ejecuta esto en un conjunto de datos de la química del suelo con las variables pH, $\mathrm{Ca^{2+}}$ TotalCarbon, podría encontrar que el primer componente era

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

y el segundo componente

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

Estos componentes se pueden seleccionar libremente entre las variables medidas, y se eligen aquellos 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 demás.

En una ordenación restringida, tenemos dos conjuntos de datos, pero no somos libres de seleccionar cualquier combinación lineal del primer conjunto de datos (los datos químicos del suelo anteriores) que queramos. En cambio, tenemos que seleccionar las combinaciones lineales de las variables del segundo conjunto de datos que mejor expliquen la variación del primero. Además, en el caso del ACP, el conjunto de datos es la matriz de respuesta y no hay predictores (se puede 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 has explicado qué variables son la respuesta, normalmente se desea explicar la variación en las abundancias o la composición de esas especies (es decir, las respuestas) utilizando las variables explicativas ambientales.

La versión restringida del ACP es una cosa que en los círculos ecológicos se llama Análisis de Redundancia (RDA). Esto supone un modelo de respuesta lineal subyacente para las especies, lo cual no es apropiado o sólo es apropiado si se tienen gradientes cortos a lo largo de los cuales las especies responden.

Una alternativa al ACP es algo llamado análisis de correspondencia (AC). No tiene restricciones, pero tiene un modelo de respuesta unimodal subyacente, que es algo más realista en términos de cómo responden las especies a lo largo de gradientes más largos. Tenga en cuenta también que los modelos de AC las abundancias relativas o la composición El PCA modela las abundancias brutas.

Existe una versión restringida de CA, conocida como restringido o análisis de correspondencia 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 el RDA como en el CCA, el objetivo es modelar la variación de la abundancia o la composición de las especies como una serie de combinaciones lineales de las variables explicativas.

Por la descripción, parece que su mujer quiere explicar la variación en la composición de las especies de milpiés (o su abundancia) en función de las otras variables medidas.

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

• a medida que se aumenta el número de variables explicativas, las restricciones son cada vez menores y no se están extrayendo componentes/ejes que expliquen la composición de las especies de forma óptima, y
• Con el CCA, a medida que se aumenta el número de factores explicativos, se corre el riesgo de inducir un artefacto de curva en la configuración de los puntos del gráfico del CCA.
• La teoría en la que se basan el RDA y el CCA está menos desarrollada que los métodos estadísticos más formales. Sólo podemos elegir razonablemente qué variables explicativas mantener utilizando 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 la regresión) y tenemos que utilizar pruebas de permutación para hacerlo.

así que mi consejo es el mismo que con la regresión; piense de antemano cuáles son sus hipótesis e incluya variables que reflejen esas hipótesis. No lo hagas sólo hay que echar todas las variables explicativas en la mezcla.

Ejemplo

Ordenación sin restricciones

PCA

Mostraré un ejemplo en el que se comparan PCA, CA y CCA utilizando el vegano paquete 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)

vegano no estandariza la inercia, a diferencia de Canoco, por lo que la varianza total es de 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 el ~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

A partir de las puntuaciones de las muestras y de las especies en los dos primeros componentes principales se puede trazar un biplot

> plot(pcfit)

Aquí hay dos cuestiones

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

CA

Una AC podría ayudar con estos dos puntos, ya que maneja mejor los gradientes largos debido al modelo de respuesta unimodal, y modela la composición relativa de las especies, no las abundancias brutas.

El vegano / R El código 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 alrededor del 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

El gráfico conjunto de las puntuaciones de las especies y del lugar está ahora menos dominado por unas pocas especies

> plot(cafit)

La elección de PCA o CA debe estar determinada por las preguntas que desee hacer a los datos. Normalmente, con los datos de las especies nos interesa más la diferencia en el conjunto de especies, por lo que el AC es una opción popular. Si tenemos un conjunto de datos de variables ambientales, por ejemplo el agua o la química del suelo, no esperaríamos que respondieran de forma unimodal a lo largo de los gradientes, por lo que el AC sería inapropiado y el ACP (de una matriz de correlación, utilizando scale = TRUE en el rda() ) sería más apropiado.

Ordenación restringida; CCA

Ahora bien, si tenemos un segundo conjunto de datos que deseamos utilizar para explicar patrones en el conjunto de datos de la primera especie, debemos utilizar una ordenación restringida. A menudo la elección es CCA, pero RDA es una alternativa, como lo es RDA después de la transformación de los datos para permitirle manejar mejor los datos de las especies.

data(varechem)                          # load explanatory example data

Reutilizamos el cca() pero suministramos dos marcos de datos ( X para las especies, y Y para las variables explicativas/predictivas) o una fórmula del modelo que enumera la forma del modelo que deseamos ajustar.

Para incluir todas las variables podríamos utilizar varechem ~ ., data = varechem como la fórmula para incluir todas las variables - pero como dije arriba, esto 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)

Por supuesto, ahora la tarea consiste en averiguar cuál de esas variables es realmente importante. Uno de los problemas de utilizar todas las variables en esta ordenación es que hemos creado una configuración arqueada en las puntuaciones de las muestras y de las especies, que es un mero artefacto de utilizar demasiadas variables correlacionadas.

Si quiere saber más sobre esto, consulte el vegano documentación o un buen libro sobre análisis de datos ecológicos multivariantes.

Relación con la regresión

Es más sencillo ilustrar el vínculo con la RDA, pero el CCA es igual, excepto que todo implica sumas marginales de filas y columnas de tablas de dos vías como pesos.

En esencia, la RDA equivale a la aplicación del ACP 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 conseguir que las puntuaciones de los ejes (cargas) coincidan, pero invariablemente éstas están escaladas (o no), así que tengo que investigar exactamente cómo se están haciendo aquí.

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