Análisis de datos ecológicos (I)

Práctica 2, Ecología II

Grado en Biología

Introducción y objetivos

En esta práctica usaremos R para realizar varios tipos de análisis de datos de uso habitual en la investigación ecológica. Utilizando diversos ejemplos de datos, tanto simulados como reales, nos centraremos en el proceso de elección de las técnicas de análisis en función de las características de los datos, en el procedimiento de diseño de los modelos y en la interpretación de los resultados (coeficientes, estadísticos, representaciones gráficas). Para profundizar más en el análisis estadístico de datos ecológicos (o biológicos en general) puedes consultar los libros de Gotelli & Ellison (2018) y Venables & Ripley (2002). También puedes realizar el seminario disponible en: https://t.um.es/aedbr. Ahí encontrarás más información, ejercicios, ejemplos y métodos de análisis más avanzados.

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( "https://webs.um.es/jfcalvo/eco2.RData" ) )

Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "aves" ) o info( "islas" ). Teclea info() para más información.

Para realizar los ejercicios propuestos se necesita la instalación adicional del paquete de R effects. Si no lo tenemos instalado utilizaremos:

install.packages( "effects" )

Una vez instalado lo cargaremos en memoria:

library( effects )

ANOVA y regresión lineal

Comenzaremos realizando análisis de la varianza (ANOVA) utilizando los datos simulados del objeto ejemplo, que contiene información sobre una variable de respuesta ecológica (cobertura de una especie de planta), una variable ambiental cualitativa (tipo de hábitat) y una variable ambiental cuantitativa (humedad del suelo).

ejemplo
str( ejemplo )

Conviene siempre hacer inicialmente una representación gráfica de los datos:

plot( cob ~ hum, data = ejemplo, col = hab, pch = 19, cex = 1.5 )

Para analizar el efecto de la variable cualitativa realizaremos un análisis de la varianza (ANOVA) con la función aov (y summary, para mostrar los resultados con más detalle):

aov( cob ~ hab, data = ejemplo ) -> m1
summary( m1 )

Para representar el resultado del modelo utilizaremos la función predictorEffects del paquete effects:

plot( predictorEffects( m1 ) )

Con la variable independiente cuantitativa haremos una regresión (modelo lineal, con errores de distribución normal o gaussiana):

lm( cob ~ hum, data = ejemplo ) -> m2
summary( m2 )
plot( cob ~ hum, data = ejemplo )
abline( 14.2693, 0.7956 )

Una representación adecuada de modelo debe incluir el intervalo de confianza. Para ello utilizaremos la función predictorEffects del paqueteeffects:

plot( predictorEffects( m2 ) )

La cobertura estimada por este modelo para un valor determinado de la variable independiente puede obtenerse aplicando la ecuación estimada. Por ejemplo, para un valor de humedad igual a 20:

14.2693 + 0.7956 * 20

También podemos usar un modelo de regresión para hacer un análisis equivalente al ANOVA:

lm( cob ~ hab, data = ejemplo) -> m3
summary( m3 )
plot( predictorEffects( m3 ) )

Ejercicio 1

El objeto pinos contiene una tabla de datos con las estimas medias de dbh de los ejemplares de pino carrasco presentes en 83 parcelas ubicadas en diferentes zonas de un gradiente de ladera (variable cualitativa zona) y con diferentes pendientes del terreno (variable cuantitativa pendiente).

• Realiza una análisis de la varianza con zona e interpreta los resultados.

summary( aov( dbh ~ zona, data = pinos ) -> mp1 )
plot( predictorEffects( mp1 ) )

• Las comparaciones múltiples (entre zonas) las podemos hacer con el test de Tukey:

TukeyHSD( mp1 )
plot( TukeyHSD( mp1 ) ) 

• Realiza un análisis de regresión lineal con pendiente e interpreta los resultados en comparación con los del ANOVA anterior.

summary( lm( dbh ~ pendiente, data = pinos ) -> mp2 )
plot( predictorEffects( mp2 ) )

Análisis de relaciones especies-área con regresión lineal

Aplicaremos un modelo de regresión lineal para analizar relaciones especies-área con los datos de riqueza de plantas en las islas Galápagos. Los datos están incluidos en el objeto islas:

str( islas )
islas[[ 1 ]]
lm( log( riqueza ) ~ log( área ), data = islas[[ 1 ]] ) -> mrsa
summary( mrsa )

Podemos hacer la representación gráfica en escala no logarítmica con:

curve( exp( 3.25520 ) * x ^ 0.3478, 1, 2500 )

También podemos usar la función predictorEffects, aunque en este caso, la representación de los intervalos de confianza requiere otro argumento:

plot( predictorEffects( mrsa ), axes = list( y = list( transform = exp, lab = "riqueza" ) ) )

Ejercicio 2

Realiza los análisis de regresión lineal, y representa las curvas de especies-área, de los datos de riqueza de aves de los archipiélagos de las Indias Occidentales (islas[[ 2 ]]) y de las Indias Orientales (islas[[ 3 ]]).

Modelos lineales generalizados

Los modelos de regresión lineal no son adecuados para determinados tipos de datos, como los conteos o los datos de presencia/ausencia (1/0). En estos casos es necesario utilizar modelos lineales generalizados (GLM). En el siguiente ejercicio analizaremos los datos de conteos del mito (Aegithalos caudatus) en relación con la variable altitud. Estos datos están incluidos en el objeto aves, una lista que contiene una matriz de datos con 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 variables ambientales de los transectos (aves$amb), y una tercera con información adicional sobre la localización de los transectos (aves$loc).

str( aves )

Generalmente las matrices de datos en R se denominan tablas de datos (data frame), reservándose el término matriz (matrix) a objetos con propiedades que permiten realizar operaciones matriciales (en el sentido matemático de la expresión). A efectos prácticos, la principal diferencia es que las tablas de datos pueden incorporar variables de diferentes características (por ejemplo, de texto y numéricas).

Para desarrollar los ejercicios más cómodamente crearemos una nueva tabla de datos combinando la matriz de conteos de aves y la de datos ambientales:

data.frame( aves$sp, aves$amb ) -> avesglm
head( avesglm )
plot( Aegithalos.caudatus ~ altitud, data = avesglm )

Usaremos la función glm con el argumento family = poisson para obtener un modelo de regresión discreta o de Poisson:

glm( Aegithalos.caudatus ~ altitud, family = poisson, data = avesglm ) -> mg1
summary( mg1 )
plot( predictorEffects( mg1 ), type = "response" )

La expresión de un modelo de regresión discreta es:

\[\ln\left( Ey \right) = b_{0} + b_{1}x \quad \rightarrow \quad Ey = e^{b_{0} + b_{1}x}\]

Esto implica que, por ejemplo, la abundancia estimada de mitos en un transecto a 1200 m s. n. m. se calcula como:

exp( 0.3339066 + 0.0005314 * 1200 )

Si queremos probar el ajuste a un modelo unimodal, solo tenemos que añadir en la ecuación el término cuadrático de la variable cuantitativa:

glm( Aegithalos.caudatus ~ altitud + I( altitud ^ 2 ), family = poisson, data = avesglm ) -> mg2
summary( mg2 )
plot( predictorEffects( mg2 ), type = "response" )

Como ejemplo de datos de presencia/ausencia seguiremos con el mito, transformando los datos de conteos en unos y ceros, y creando una nueva tabla de datos:

data.frame( ( aves$sp > 0 ) * 1, aves$amb ) -> avespres
head( avespres )
plot( Aegithalos.caudatus ~ altitud, data = avespres )

Usaremos de nuevo la función glm, pero ahora con el argumento family = binomial, para obtener un modelo de regresión logística:

glm( Aegithalos.caudatus ~ altitud, family = binomial, data = avespres ) -> mg3
summary( mg3 )
plot( predictorEffects( mg3 ), type = "response", ylim = c( 0, 1 ) )

Este tipo de modelos estiman probabilidades de presencia (u ocupación) mediante el ajuste a una curva logística, cuya expresión es:

\[\ln\left( \frac{p}{1 - p} \right) = b_{0} + b_{1}x \quad \rightarrow \quad p = \frac{e^{b_{0} + b_{1}x}}{1 + e^{b_{0} + b_{1}x}}\]

La obtención de la probabilidad de presencia estimada para un determinado valor de la variable explicativa (altitud) se obtiene fácilmente con la función plogis. Por ejemplo:

plogis( -1.4242457 + 0.0014181 * 1200 )

Ejercicio 3

Prueba el ajuste de los datos de presencia/ausencia del mito a una curva logística unimodal:

glm( Aegithalos.caudatus ~ altitud + I( altitud ^ 2 ), family = binomial, data = avespres) -> mg4
summary( mg4 )
plot( predictorEffects( mg4 ), type = "response", ylim = c( 0, 1 ) )

Evaluación

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

Bibliografía

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

• Venables WN, Ripley BD (2002) Modern Applied Statistics with S. 4ª ed. Springer, New York.

Descripción de los datos

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