Para esta lección volveremos a la base de datos de abundancia de especies. Recordemos que para importar la base de datos necesitamos conocer la ruta donde se ubica el archivo y su nombre. En mi caso:

spp <- read.csv("Datos-bases/Datos-sp.csv")

El primer examen que haremos es una gráfica de dispersión entre pares de variables, para poder identificar cuáles están correlacionadas entre sí, positiva o negativamente. Usaremos la función pairs, la cual produriá automáticamente el gráfico para cada par de variables (figura ??):

pairs(spp, upper.panel = NULL)

Este análisis gráfico podemos acompañarlo de una tabla con los coeficientes de correlación estimados (tabla ??):

cor(spp)
##            Sp1         Sp2         Sp3         Sp4        Sp5
## Sp1  1.0000000  0.83521375 -0.06045670 -0.11719348  0.1971041
## Sp2  0.8352138  1.00000000 -0.04319937 -0.05909854  0.1573974
## Sp3 -0.0604567 -0.04319937  1.00000000  0.78338539 -0.7644125
## Sp4 -0.1171935 -0.05909854  0.78338539  1.00000000 -0.8196405
## Sp5  0.1971041  0.15739739 -0.76441254 -0.81964054  1.0000000

Para representar la relación entre tres variables simultáneamente necesitamos un gráfico en 3d:

library(lattice)
cloud(Sp3 ~ Sp4 + Sp5, spp, type="p",screen=c(x=-55,y=-30,z=-20),zlab=list(rot=90), pch=16)
Gráfico de dispersión para tres variables especies correlacionadas entre sí.

Figure 0.1: Gráfico de dispersión para tres variables especies correlacionadas entre sí.

En este caso decidí mostrar el gráfico de dispersión entre Sp3, 4 y 5 porque están mutuamente correlacionadas (figura ??). Por lo tanto, podríamos representar a estas tres especies como una sola en una comunidad, pues se predicen mutuamente. No obstante, los grupos de variables que se predicen mutuamente no siempre se pueden identificar tan claramente como en este caso. Los análisis multivariados nos ayudarán a identificar grupos de variables que se pueden resumir para disminuir la complejidad de una base de datos.

1 Análisis de componentes principales

Una técnica comunmente utilizada es el análisis de componentes principales, el cual por medio de la correlación entre variables genera componentes que representan la mayor proporción de variabilidad en el menor número de variables. En R, podemos hacer un análisis de componentes principales con la función prcomp:

pca <- prcomp(spp)

El objeto resultante de este análisis podemos graficarlo, para ver cuánta varianza contiene cada uno de los componentes que genera y cómo están constituidos:

plot(pca)
Gráfico de la varianza que contiene cada componente.

Figure 1.1: Gráfico de la varianza que contiene cada componente.

Este gráfico contiene en el eje de las \(x\) la identificación del componente (1-5), y en el \(y\) la varianza. Los componentes simbre están ordenados descendentemente con base en la varianza. Generalmente se recomienda utilizar los componentes antes del “codo” que se forma entre los dos ejes. En este caso resulta evidente que los componentes más relevanes son el 1 y 2, lo cual podemos revisar con el porcentaje de varianza que contienen:

pca$sdev^2/sum(pca$sdev^2)
## [1] 0.52810535 0.35714445 0.04934741 0.03490572 0.03049707

Ahora podemos ver cómo están constituidos estos dos componentes con el gráfico:

biplot(pca)

Para interpretar este gráfico, podemos guiarnos con la figura ??. En ?? vemos que las flechas rojas con las etiquetas Sp1 y Sp2 apuntan en una dirección muy similar, lo cual quiere decir que están positivamente correlacionadas, y lo mismo ocurre con Sp3 y Sp4. Estos dos grupos de especies, Sp1-Sp2 y Sp3-Sp4 son ortogonales puesto que no están correlacionados, lo cual podemos ver con las flechas prácticamente perpendiculares entre sí. En términos prácticos podemos inferir que la presencia de Sp1 y Sp2 no afecta la presencia de Sp3 y Sp4. En relación a Sp5, su vector apunta en dirección opuesta de los vectores de Sp3 y Sp4, debido a que están negativamente correlacionadas (Figura ??).

2 Mediciones de similitud

Así como hemos analizado en mayor detalle la relación multivariada entre todas las especies que conforman la base spp, es posible compara entre sí todos los sitios de muestreo, de modo que podamos identificar cuánto se parecen entre sí. Este tipo de comparaciones se tiene que hacer forzosamente entre pares de sitios o unidades de muestreo, y por lo tanto, seleccionaremos aleatoriamente 10 sitios de los 100 contenidos en spp:

set.seed(10)
spp.10 <- spp[sample(1:100, 10),]

knitr::kable(spp.10, caption = "Los 10 sitios seleccionados al azar para calcular los índices de similitud.")
Table 2.1: Los 10 sitios seleccionados al azar para calcular los índices de similitud.
Sp1 Sp2 Sp3 Sp4 Sp5
9 6 10 20 11 4
74 4 10 20 12 4
76 7 11 19 8 6
55 5 10 19 9 7
72 6 11 20 10 6
54 6 11 18 9 6
39 6 11 19 9 5
83 5 10 19 10 6
88 5 11 18 9 6
15 5 9 20 10 4

Uno de los índices más usados en ciencias para medir la similitud entre sitios es el de Jaccard. Su implementación es muy sencilla por lo que hay pocas implementaciones en paquetes de R. En índice de similitud entre dos sitios es:

\[J(A, B) = \frac{|A \cap B|}{|A \cup B|}\]

donde \(|A \cap B|\) es el número de elementos de \(A\) que también se encuentran en \(B\), y \(|A \cup B|\) es el número de elementos que están en ambos conjuntos.

Escribiremos una función para calcular el índice de Jaccard para un par de sitios cualquiera:

jaccard <- function(A, B){
    interseccion = length(intersect(A, B))
    union = length(A) + length(B) - interseccion
    return (interseccion/union)
}

Entonces para calcular el índice para los sitios 1 y 2:

jaccard(spp.10[1,], spp.10[2,])
## [1] 0.4285714

Esta operación debemos repetirla para todos los pares de sitios, cuyo resultado se muestra en la figura 2.1

m <- matrix(0, 10, 10)
for (i in 1:10) {
    for(j in 10:1){
        m[i, j] <- jaccard(spp.10[i, ], spp.10[j, ])
    }
}
image(m)
Índice de jaccard para todos los pares de sitios de `spp.10`.

Figure 2.1: Índice de jaccard para todos los pares de sitios de spp.10.

2.1 Mediciones más generales de similitud

La similitud medida con el índice de Jaccard está basada en la repetición de elementos en dos conjuntos, y debe ser derivada para todos los pares de sitios. En posible, sin embargo, encontrar medidas de similitud en relación a tendencias inherentes a la base de datos que se está analizando. Una de estas medidas es la distancia. Así como podemos medir distancias físicas entre dos puntos, es posible medir distancias entre conjuntos de variables. Como vimos en la sección de PCA, los datos pueden contener cierta estructura, que se manifiesta como correlación. Paa medir distancias es necesario tomar en cuenta dicha correlación. Para medir distancia entre las características de un sitio y las tedencias encontradas en toda la base de datos pueden hacerse en relación al promedio aritmético de cada una de las variables medidas. Entonces, en la base de datos spp.10 podemos calcular la media de cada una de las columnas con la función colMeans:

medias <- colMeans(spp.10)

La distancia entre cada sitio de muestreo y el vector de medias puede medirse con el teorema de Pitágoras, lo cual se llama distancia Euclidiana. En este caso, para tomar en cuenta la correlación entre Sp1-Sp2, y Sp3-Sp4, necesitamos una matriz de covarianza:

mat.cov <- cov(spp.10)
mat.cov
##            Sp1        Sp2        Sp3        Sp4        Sp5
## Sp1  0.7222222  0.3333333 -0.1111111 -0.6111111  0.2222222
## Sp2  0.3333333  0.4888889 -0.3111111 -0.4222222  0.3777778
## Sp3 -0.1111111 -0.3111111  0.6222222  0.6222222 -0.5333333
## Sp4 -0.6111111 -0.4222222  0.6222222  1.3444444 -0.8666667
## Sp5  0.2222222  0.3777778 -0.5333333 -0.8666667  1.1555556

Una vez que tenemos la matriz de covarianza, calculamos la distancia entre las abundancias de cada sitio y la media de cada una de las abundancias de cada especie:

dists <- mahalanobis(spp.10, center = medias, cov = mat.cov)
dists.df <- data.frame(Sitio = paste0("Sitio ", 1:10),
                       Distancia = dists)

knitr::kable(dists.df, caption = "Distancias de Mahalanobis al centroide de las abundancias registradas.")
Table 2.2: Distancias de Mahalanobis al centroide de las abundancias registradas.
Sitio Distancia
9 Sitio 1 5.275190
74 Sitio 2 5.540609
76 Sitio 3 3.747981
55 Sitio 4 4.813868
72 Sitio 5 5.538502
54 Sitio 6 4.215623
39 Sitio 7 3.212931
83 Sitio 8 2.399824
88 Sitio 9 3.937566
15 Sitio 10 6.317905

3 Mediciones de diversidad

Al medir qué tan equitativamente están distribuidos los recursos entre las especies que conforman un ensable, es necesario derivar algún número que permita la comparación entre comunidades, tomando en cuenta tanto las especies presentes como sus abundancias. Aquí, entran lo que se conoce como índices de diversidad, y son medidas multivariadas de equitatividad. Dos de los índices más utilizados son el de Simpson y el Shanon-Wiener. Aquí calcularemos el índice de Simpson para cada uno de los 10 sitios de muestreo. El paquete abdiv tiene una implementación de dicho índice:

library(abdiv)
## 
## Attaching package: 'abdiv'
## The following object is masked _by_ '.GlobalEnv':
## 
##     jaccard
simp <- apply(spp.10, 1, simpson)
simp.df <- data.frame(Sitio = 1:10, Simp.ind = simp)

knitr::kable(simp.df, caption = "Índice de diversidad de Simpson calculado para los sitios 1 a 10.")
Table 3.1: Índice de diversidad de Simpson calculado para los sitios 1 a 10.
Sitio Simp.ind
9 1 0.7412534
74 2 0.7296000
76 3 0.7574010
55 4 0.7536000
72 5 0.7532930
54 6 0.7608000
39 7 0.7504000
83 8 0.7512000
88 9 0.7555185
15 10 0.7300347

Código fuente de este tutorial

Regresar al índice del curso