load(url( "https://webs.um.es/jfcalvo/eco2.RData" ) )
Análisis de datos ecológicos (I)
Práctica 2, Ecología II, Grado en Biología
Actualizado el 17/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 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.
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:
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 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 )(B) ANOVA, regresión lineal y ANCOVA
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( ejemplo$cob ~ ejemplo$hum, col = ejemplo$hab )O mejor:
plot( cob ~ hum, data = ejemplo, col = hab, pch=19 )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 * 20Los modelos lineales nos permiten además combinar variables de diferente naturaleza. Por ejemplo, un análisis de la covarianza (ANCOVA) proporciona el siguiente resultado:
lm( cob ~ hab + hum, data = ejemplo ) -> m3
summary( m3 )
plot( cob ~ hum, data = ejemplo, col = hab, pch = 19 )
abline( 23.82130, 0.75231 )
abline( 23.82130 -17.24897, 0.75231, col = 2 )Finalmente interpretaremos un modelo más complejo que incluye términos de interacción entre las variables:
lm( cob ~ hab * hum, data = ejemplo ) -> m4
summary( m4 )
plot( cob ~ hum, data = ejemplo, col = hab, pch = 19 )
abline( 31.71022, 0.39647 )
abline( 31.71022 -30.84935, 0.39647 +0.63256, col = 2 )La representación con los intervalos de confianza la haríamos con:
plot( predictorEffects( m4, ~ hum ), multiline = TRUE, confint = list( style = "bands" ) )(C) 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:
names( islas )
islas[[ 1 ]]
lm( log( riqueza ) ~ log( área ), data = islas[[ 1 ]] ) -> m5
summary( m5 )Podemos hacer la representación gráfica en escala no logarítmica con:
curve( exp( 3.25520 ) * x ^ 0.3478, 1, 2500 )En este caso, la representación de los intervalos de confianza requiere otro argumento:
plot( predictorEffects( m5 ), axes = list( y = list( transform = exp ) ) )(D) 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( avesglm[ , 1 ] ~ 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( avesglm[ , 1 ] ~ 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( avesglm[ , 1 ] ~ altitud + I( altitud ^ 2 ), family = poisson, data = avesglm ) -> mg2
summary( mg2 )
plot( predictorEffects( mg2, xlevels = 100 ), type = "response" )El valor proporcionado con argumento xlevels permite obtener una representación “suavizada” de la curva.
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( avespres[ , 1 ] ~ 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( avespres[ , 1 ] ~ 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 )También podemos probar el ajuste a curvas logísticas unimodales, aunque en este caso no es estadísticamente significativo:
glm( avespres[ , 1 ] ~ altitud + I( altitud ^ 2 ), family = binomial, data = avespres) -> mg4
summary( mg4 )
plot( predictorEffects( mg4, xlevels = 100 ), type = "response", ylim = c( 0, 1 ) )(E) Tablas de contingencia
Las tablas de contingencia constituyen la forma más habitual de presentación de datos en los que tanto la variable de respuesta como la variable o variables explicativas son de naturaleza cualitativa. Usaremos como ejemplo los datos procedentes de un muestreo de presencia/ausencia de carrasca (Quercus rotundifolia) en cuadrículas UTM de 250 x 250 m realizado en la Sierra de la Torrecilla (Región de Murcia), en el que cada cuadrícula fue asignada a una determinada categoría de uso o tipo de hábitat (agrícola, forestal o mixto).
carrascaEl análisis de la tabla puede realizarse mediante pruebas clásicas, como el test de \(\chi^2\), que permite evaluar la independencia entre ambas variables:
chisq.test( carrasca )No obstante, es más recomendable elaborar las tablas de contingencia con formato de tablas de frecuencias. La conversión la podemos hacer con la función freq.table, incluida en el archivo eco2.RData:
freq.table( carrasca, "uso", "estado" ) -> carrascaF
carrascaFCon este formato pueden aplicarse diversos tipos de modelos de regresión que proporcionan mayor flexibilidad para los análisis, incluidos los de tablas de contingencia con más de un factor. Por ejemplo, aplicaremos a la nueva tabla un análisis de regresión logística y un análisis loglineal (con regresión de Poisson), que proporcionan resultados equivalentes:
glm( estado ~ uso, weights = Freq, family = binomial, data = carrascaF ) -> mf1
summary( mf1 )
glm( Freq ~ uso * estado, family = poisson, data = carrascaF ) -> mf2
summary( mf2 )En este último modelo, los únicos tests de interés son los que corresponden a los coeficientes de los términos de interacción, que determinan la significación estadística de diferencias entre el uso de referencia (intercept) y el resto. El test global para ambos análisis se realiza con la función anova (no confundir con aov):
anova( mf1, test = "Chisq" )
anova( mf2, test = "Chisq" )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.
Venables WN, Ripley BD (2002) Modern Applied Statistics with S. 4ª ed. Springer, New York.
Descripción de los datos y funciones
Utiliza la función info. Por ejemplo: info("ejemplo") o info("freq.table").