Análisis de datos ecológicos (II)

Práctica 3, Ecología II, Grado en Biología

Actualizado el 26/10/2023

Profesor: José Francisco Calvo Sendín
Área de Ecología, Facultad de Biología, Universidad de Murcia
jfcalvo@um.es | https:/webs.um.es/jfcalvo

1. Introducción y objetivos

En esta práctica realizaremos varios ejercicios, a través del uso de una matriz de datos reales, para aprender diversos procedimientos relacionados con el manejo de tablas de datos ecológicos en R. Realizaremos también cálculos y análisis descriptivos de los datos, así como representaciones gráficas de tipo comparativo entre variables, de relaciones entre variables biológicas y ambientales, y de carácter espacial. Finalmente aplicaremos técnicas de análisis multivariante (ordenación y clasificación), y realizaremos un ejemplo de cálculo de diversidad funcional. Puedes encontrar más información sobre la aplicación de métodos multivariantes a datos ecológicos en el libro de Gotelli y Ellison (2018).

2. Desarrollo de la sesión

(A) Preparación

Los datos necesarios para el desarrollo de la práctica se encuentran disponibles en el archivo eco2.RData. Una vez iniciado R cargaremos este archivo desde el servidor de la asignatura:

load( url( "http://webs.um.es/jfcalvo/eco2.RData" ) )

Usando ls() podemos ver los datos cargados. Utiliza la función info para más información. Por ejemplo: info() o info("aves").

Para realizar los ejercicios propuestos se necesitan dos paquetes de R que requieren una instalación adicional (terra y vegan). Si no los tenemos instalados utilizaremos:

install.packages( c( "terra", "vegan") )

Una vez instalados los cargaremos en memoria:

library( terra )
library( vegan )

(B) Manejo de tablas de datos de comunidades

Trabajaremos con las tablas de datos del objeto aves, que conocemos de la práctica 2. Recordemos que se trata de una lista que contiene una matriz o tabla de datos de conteos de individuos de 16 especies de aves (aves$sp) en 140 transectos (itinerarios de censo de 1 km de longitud) en la Región de Murcia, otra con datos de variables ambientales de los transectos (aves$amb), y otra con información adicional sobre la localización de los transectos (aves$loc):

str(aves)

Para realizar determinados ejercicios resulta cómodo unir las tres matrices en una única tabla de datos, usando la función data.frame (o alternativamente cbind):

data.frame( aves$sp, aves$amb, aves$loc ) -> aves2
head( aves2 )

A menudo es necesario seleccionar partes de la matriz, en función, por ejemplo, de alguna de las variables ambientales. Así, si quisiéramos trabajar solo con los transectos de las sierras litorales, usaríamos:

subset( aves2, aves2$litoral == 1 ) -> aveslitoral
aveslitoral

Si necesitáramos usar más de una condición para la selección, emplearíamos operadores lógicos. Por ejemplo:

subset( aves2, aves2$bosque == 100 & aves2$temperatura > 15 )
subset( aves2, aves2$precipitación > 450 | aves2$sierra == "Espuña" )
subset( aves2, aves2$litoral != 0 | aves2$sierra == "Carche" )
Ejercicio 1

Selecciona los transectos del objeto aves2 localizados a altitudes menores de 500 m s. n. m. y con superficies de bosque por debajo del 50 %. ¿Qué dimensiones tiene esta nueva tabla de datos?

La transformación de variables es una tarea habitual en la gestión de matrices de datos. Por ejemplo, si nos interesara identificar los transectos con una elevada proporción de bosque (≥ 75 %), podemos crear fácilmente una nueva variable cualitativa y añadirla a la tabla de datos:

aves2$bosque >= 75
aves2$bosque >= 75 -> aves2$bosque2
head( aves2 )

También se pueden realizar codificaciones más complejas, aunque el procedimiento es más laborioso y requiere unas sencillas líneas de programación en R (ejercicio 2).

Ejercicio 2

Comprueba el efecto de las siguientes órdenes para una nueva codificación de la variable bosque:

aves2$bosque3 <- NA
for ( i in 1:140 ) if ( aves2$bosque[ i ] <= 25 ) aves2$bosque3[ i ] <- "baja"
for ( i in 1:140 ) if ( aves2$bosque[ i ] > 25 & aves2$bosque[ i ] < 75 ) aves2$bosque3[ i ] <- "media"
for ( i in 1:140 ) if ( aves2$bosque[ i ] >= 75 ) aves2$bosque3[ i ] <- "alta"

(C) Cálculos

Además de las funciones de sumas de filas (rowSums) y columnas (colSums), R dispone de otras funciones para realizar diversos tipos de cálculos sobre los datos de las tablas. La función apply permite aplicar varias funciones sobre filas y columnas. Con esta función hay que indicar con un 1 si queremos trabajar con las filas, o con un 2 si queremos trabajar con las columnas. Por ejemplo:

apply( aves$sp, 1, sum )
apply( aves$sp > 0, 1, sum ) -> aves2$riqueza
head( aves2 )
apply( aves$sp, 2, mean )
apply( aves$sp, 2, var )
apply( aves$sp, 2, table )
Ejercicio 3

Selecciona los transectos del objeto aves2 localizados a altitudes mayores de 1000 m s. n. m. y con temperaturas medias anuales por debajo de 14 ºC.

  1. ¿Cuántos individuos de las distintas especies de aves se observaron en total en los transectos seleccionados?
  2. ¿Cuál es el número medio de individuos de las especies del género Curruca en esos transectos?
  3. ¿En qué transecto se registró la mayor riqueza de especies?

Una función muy útil es aggregate, que permite realizar cálculos sobre las variables en función de uno o varios criterios de agrupación. Por ejemplo, la suma total de individuos de las 16 especies de aves contadas en los itinerarios de cada sierra se puede calcular con:

aggregate( aves$sp, list( aves$loc$sierra ), sum )

Para conocer el número máximo de individuos contados en itinerarios de sierras litorales y no litorales, en altitudes por encima y por debajo de los 300 m s. n. m., usaríamos:

aggregate(aves$sp, list( aves2$litoral, aves2$altitud > 300 ), max )
Ejercicio 4

Calcula el número total de individuos de cada especie contados en itinerarios con porcentajes de superficie de bosque por encima y por debajo del 75 %.

(D) Representaciones gráficas

Las representaciones gráficas aportan información muy relevante sobre las características de los datos. Para datos de conteos son muy útiles los diagramas de barras que representan las frecuencias de números de individuos registrados en la muestra. Usaremos la función barplot, aplicada sobre una tabla de frecuencias. Por ejemplo, para Parus major:

table( aves2$Parus.major )
barplot( table( aves2$Parus.major ) )

Más habituales son las representaciones de datos biológicos en función de variables ambientales. Por ejemplo, podemos examinar la distribución de la abundancia de Serinus serinus, en función de la altitud:

plot( Serinus.serinus ~ altitud, data = aves2 )

La representación puede realizarse utilizando distintos tipos de símbolos y diferentes colores, en función, por ejemplo, de los valores de una variable cualitativa:

plot( Serinus.serinus ~ altitud, data = aves2, pch = 19, col = litoral )

También podemos representar conjuntamente datos de dos o más especies, e incluso añadir curvas polinómicas ajustadas a los datos. Por ejemplo:

points( Fringilla.coelebs ~ altitud, data = aves2, pch = 19, col =4 )
lines( smooth.spline( aves2$Serinus.serinus ~ aves2$altitud ), lwd = 2 )
lines( smooth.spline( aves2$Fringilla.coelebs ~ aves2$altitud), lwd = 2, col =4 )

Por otra parte, la comparación de las abundancias medias de las distintas especies entre sí puede realizarse utilizando la función boxplot, que representa diagramas de caja (box plots):

boxplot( aves$sp )

En un box plot se representan la mediana, el rango intercuartílico, el rango (máximo y mínimo) y los valores extremos (outliers) de las variables.

La función boxplot también se puede usar en combinación con variables ambientales cualitativas:

boxplot( Serinus.serinus ~ litoral, data = aves2 )
boxplot( Serinus.serinus ~ litoral * (bosque > 50), data = aves2 )

En el ejemplo de las aves, al disponer de información espacial (coordenadas geográficas de los itinerarios), también podemos realizar representaciones gráficas en forma de mapas. Usaremos varias funciones del paquete terra.

Primero representaremos el mapa con los límites de la Región de Murcia (objeto murcia incluido en el archivo eco2.RData), y sobre él la localización de los transectos (argumento add = TRUE):

plot( vect( murcia ) )
vect( aves2, geom = c( "UTMx", "UTMy" ) ) -> avesmap
plot( avesmap, add = TRUE, cex = 0.5 )

Para visualizar la distribución de la abundancia de una determinada especie podemos representar el tamaño del punto de cada transecto en función del número de individuos de esa especie contados en cada transecto:

plot( avesmap, add = TRUE, col = 6, cex = aves2$Curruca.undata, alpha = 0.2 )
Ejercicio 5

Representa la distribución de la abundancia de Lophophanes cristatus sobre el mapa de la Región de Murcia.

(E) Ordenación

Las técnicas de ordenación, junto con las técnicas de clasificación (sección F), constituyen un amplio grupo de métodos de análisis denominados multivariantes. A diferencia de los métodos utilizados en la práctica anterior, en lugar de una única variable de respuesta, las técnicas multivariantes analizan conjuntamente múltiples variables (habitualmente múltiples especies), por lo que resultan apropiadas para el estudio de los patrones de respuesta de las comunidades a variables ambientales.

Utilizaremos varias funciones del paquete vegan y trabajaremos con las matrices del objeto aves. Realizaremos inicialmente un análisis de correspondencias, apropiado para datos procedentes de conteos, con la función cca:

cca( aves$sp ) -> ord1
summary( ord1 )
plot( ord1, type = "text" )

Como estamos analizando muchas unidades de muestreo, la representación gráfica no permite apreciar adecuadamente la disposición de las especies en el plano definido por los ejes 1 y 2. Para seleccionar solo las especies usaremos:

plot( ord1, display = "species" )

La versión del análisis de correspondencias que permite analizar conjuntamente la matriz de especies y la matriz de variables ambientales se denomina correspondencias canónicas. Se realiza también con la función cca. Representaremos los resultados en una nueva ventana para poder comparar ambos análisis, teniendo en cuenta que la dirección de los ejes (positivo/negativo) es irrelevante:

x11()
cca( aves$sp, aves$amb) -> ord2
summary( ord2 )
plot( ord2 )

En este caso aparecen flechas que indican las relaciones de las variables ambientales con los ejes representados (este tipo de gráfica se denomina biplot). Para aclarar la representación podemos utilizar:

plot( ord2, display = c( "species", "bp" ) )

Otro tipo de análisis multivariante muy utilizado en ecología es el escalamiento multidimensional no métrico (nonmetric multidimensional scaling, NMDS), que calcula previamente una matriz de distancias o disimilitudes entre las variables. Como ejemplo usaremos el objeto rapaces, una lista en la que se proporciona información cualitativa sobre las características ecológicas (que consideraremos en la sección G como rasgos funcionales) de especies de aves rapaces de la Región de Murcia (rapaces$rasgos), e inventarios de presencia/ausencia de dichas especies en los Parques Regionales de la red de espacios protegidos (rapaces$parques). El análisis se realiza con la función metaMDS del paquete vegan:

rapaces
metaMDS( rapaces$rasgos ) -> ord3
plot( ord3, type = "text" )
Ejercicio 6
  1. Realiza un análisis de correspondencias con la tabla de datos rapaces$rasgos y compáralo con el escalamiento multidimensional no métrico realizado anteriormente.
  2. Realiza un escalamiento multidimensional no métrico con la tabla de datos rapaces$parques.
  3. Realiza un escalamiento multidimensional no métrico con la tabla de datos aves$sp y compáralo con el análisis de correspondencias realizado anteriormente.

(F) Clasificación

Las técnicas de clasificación permiten establecer relaciones jerárquicas y grupos (clusters), generalmente entre las unidades de muestreo (filas de la matriz), pero también entre las variables (columnas). Las clasificaciones se aplican sobre matrices de distancias o disimilitudes calculadas previamente. Es muy habitual también realizar clasificaciones utilizando los resultados previos (las coordenadas en los ejes) de un análisis de ordenación. Por ejemplo, en el caso del análisis de correspondencias de la matriz de pájaros (ord1) podemos utilizar las coordenadas de las especies para realizar una clasificación con la función hclust juntamente con dist (ambas incluidas por defecto en R). Usaremos solo los cuatro primeros ejes, que explicaban más del 50 % de la varianza (inercia):

hclust( dist( ord1$CA$v[ , 1:4 ] ) ) -> clus1
plot( clus1 )

Existen diversos métodos, tanto de agrupación para realizar la clasificación (con hclust), como para el cálculo de distancias (con dist), y que pueden proporcionar resultados muy diferentes. Por ejemplo:

x11()
hclust( dist( ord1$CA$v[ , 1:4 ], method = "canberra"), method = "average" ) -> clus2
plot( clus2 )

Si queremos establecer grupos concretos, debemos “cortar” el árbol con la función cutree. Podemos elegir una altura de corte (argumento h) o un número de grupos determinado (argumento k):

cutree( clus2, h = 2.5 )
cutree( clus2, k = 4 )

En el caso de datos binarios (1/0), como la matriz rapaces$rasgos, el método más apropiado de cálculo de distancias es binary (equivalente al índice de disimilitud de Jaccard):

hclust( dist( rapaces$rasgos, method = "binary" ) ) -> clus3
plot( clus3, hang = -1 )
Ejercicio 7

Realiza una clasificación con la tabla de datos rapaces$rasgos, utilizando la función betadiver del paquete vegan (seminario 1) para calcular el índice de disimilitud de Jaccard. Comprueba que el resultado es el mismo que el obtenido con la clasificación realizada anteriormente (clus3).

(G) Diversidad funcional

Realizaremos un ejemplo de cálculo de diversidad funcional aplicado sobre los inventarios de presencia de aves rapaces en los Parques Regionales de la Región de Murcia. Aprovechando la clasificación realizada en el apartado anterior, utilizando los rasgos funcionales de las especies, calcularemos para cada parque el índice FGR (functional group richness), definido como la longitud total de las ramas del dendrograma de rasgos que conecta las especies del inventario (Petchey y Gaston 2006). Se calcula con la función treedive del paquete vegan:

treedive( rapaces$parques, clus3 )

Es interesante comprobar cómo se obtienen los valores del índice a partir de las distancias del dendrograma, que pueden deducirse de la información contenida en objeto resultado de la clasificación.

clus3$merge
clus3$height

Por ejemplo, para simplificar los cálculos, podemos calcular la diversidad funcional de un nuevo inventario con solo dos especies: la 1 (águila calzada) y la 7 (azor común):

rapaces$parques[ 8, ] <- c( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 )
treedive( rapaces$parques, clus3 )

Según indica la primera fila de clus3$merge, el águila calzada (especie 1) se une a la especie 15 para formar el grupo 1. Por su parte, el azor (especie 7) se une a la especie 11 para formar el grupo 5 (fila 5 de clus3$merge), y a continuación se unen a la especie 16 para formar el grupo 7 (fila 7 de clus3$merge). Los grupos 1 y 7 se unen para formar el grupo 13, al que le corresponde una distancia de 0,6666667 (décimo tercer valor de clus3$height), que hay que multiplicar por 2, para obtener el valor del índice del nuevo inventario (1,333333).

3. Evaluación

Realiza la prueba de evaluación disponible en la herramienta Exámenes del Aula Virtual.

4. Bibliografía

  • Gotelli NJ, Ellison AM (2018) A Primer of Ecological Statistics. 2ª ed. Sinauer, Sunderland, MA.

  • Petchey OL, Gaston KJ (2006) Functional diversity: back to basics and looking forward. Ecology Letters 9: 741-758.

Descripción de los datos

Utiliza la función info. Por ejemplo: info("aves") o info("rapaces").