Análisis de datos ecológicos (I)

Práctica 2, Ecología II

Grado en Biología

Introducción

En esta práctica realizaremos varios tipos de análisis de regresión 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, donde encontrarás más información, ejemplos, ejercicios 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 )

Regresión lineal (y ANOVA)

Comenzaremos realizando análisis de regresión lineal 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 cuantitativa (humedad del suelo) y una variable ambiental cualitativa (tipo de hábitat).

ejemplo
str( ejemplo )

Siempre conviene 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 cuantitativa realizaremos un análisis de regresión (modelo lineal, con errores de distribución normal o gaussiana):

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

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

plot( predictorEffects( m1 ) )

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 analizar el efecto de la variable cualitativa mediante una regresión lineal, aunque en estos casos la interpretación de los resultados es algo diferente:

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

El resultado anterior es el mismo que el de un análisis de la varianza (ANOVA), que en R se realiza con la función aov (y summary, para mostrar los resultados con más detalle):

summary( aov( cob ~ hab, data = ejemplo ) )

Un análisis más completo de los datos de la tabla ejemplo podría incluir en el mismo modelo ambas variables explicativas y su interacción, aunque en este caso la interpretación de los coeficientes es más compleja.

lm( cob ~ hum * hab, data = ejemplo) -> m3
summary( m3 )
plot( predictorEffects( m3, ~ hum ), multiline = TRUE, confint = list( style = "bands" ), 
      main = "", lattice = list( key.args = list( x = 0.8, y = 0.2 ) ) )
points( cob ~ hum, data = ejemplo, col = hab, pch = 19, cex = 1.5 )

Puedes obtener más información sobre modelos con interacciones en https://t.um.es/aedbr.

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

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

TukeyHSD( mp1 )
plot( TukeyHSD( mp1 ) ) 

• Realiza un análisis de regresión lineal con la 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 )
points(islas[[1]]$área, islas[[1]]$riqueza)

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: regresión de Poisson

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 sobre características morfométricas de las especies (aves$rasgos).

str( aves )

Para desarrollar los ejercicios 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" )

La regresión de Poisson no es el único método apropiado para el análisis de datos de conteos. De hecho, los modelos de regresión binomial negativa suelen proporcionar mejores resultados. Puedes obtener más información en https://t.um.es/aedbr.

Modelos lineales generalizados: regresión logística

Para analizar datos de ocupación (presencia/ausencia: 1/0) recurriremos a modelos de regresión logística. Seguiremos en este caso con el ejemplo del 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 )

El ajuste de los datos de presencia/ausencia del mito a una curva logística unimodal lo realizaremos con:

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" ).