Introducción al diseño e interpretación de modelos ecológicos

Práctica 1, Métodos en Biología de la Conservación

Máster en Áreas Protegidas, Recursos Naturales y Biodiversidad

Introducción

En esta práctica recordaremos algunas nociones básicas de análisis estadístico y nos centraremos en el diseño e interpretación de modelos ecológicos, principalmente a través de la ejecución de diversas técnicas de análisis de regresión. Utilizaremos métodos de regresión apropiados para el análisis de datos ecológicos y ambientales de diferente naturaleza (cuantitativos, discretos, cualitativos). Pondremos especial énfasis en la interpretación de los modelos a través del análisis de las ecuaciones y sus coeficientes, así como de la representación gráfica de los resultados. Finalmente aprenderemos a realizar procedimientos de selección de modelos (basados en el criterio de información de Akaike, AIC) e inferencia multimodelo, aspecto fundamental para el desarrollo de las siguientes prácticas de la asignatura.

Para profundizar más en el análisis estadístico de datos ecológicos puedes consultar los libros de Gotelli & Ellison (2004) 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

Para la realización de los ejercicios de esta práctica utilizaremos R, un potente software estadístico y gráfico que puede descargarse gratuitamente desde la página https://cran.r-project.org/. A través de la web de la red CRAN (The Comprehensive R Archive Network) puede encontrarse la documentación necesaria para iniciarse y profundizar en el uso de R, así como de sus bibliotecas o paquetes de funciones, creadas por una amplísima comunidad de programadores e investigadores de todos los ámbitos de las ciencias. También puedes realizar el seminario disponible en https://t.um.es/eco1s1.

Una vez iniciado R debemos cargar el archivo de datos y funciones de la asignatura (mbc.RData) accediendo al servidor con la siguiente función:

load( url( "https://webs.um.es/jfcalvo/mbc.RData" ) )

Este archivo contiene todos los datos de ejemplo y las funciones necesarias para el desarrollo de los ejercicios propuestos en las prácticas de la asignatura. Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "ejemplo" ) o info( "territorios" ). Teclea info() para más información.

Como alternativa podemos descargar el archivo en nuestro ordenador accediendo con un navegador de internet a https://www.um.es/docencia/emc/mbc.RData. A continuación, una vez iniciado R, es conveniente cambiar el directorio de trabajo a dicha carpeta (> Archivo > Cambiar dir…), de forma que podamos acceder fácilmente al archivo descargado: > Archivo > Cargar área de trabajo… También podemos clicar sobre el archivo para que se abra directamente en el programa.

Para realizar esta práctica se necesitan dos paquetes de R que requieren una instalación adicional (effects y MuMIn). Si no los tenemos instalados usaremos:

install.packages( c( "effects", "MuMIn" ) )

Una vez instalados los cargaremos en memoria:

library( effects )
library( MuMIn )

ANOVA, regresión lineal y ANCOVA

Comenzaremos realizando análisis de la varianza (ANOVA) utilizando los datos simulados contenidos en el objeto ejemplo:

ejemplo

Podemos conocer la estructura de este objeto y las características de sus variables con:

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 )
model.tables( m1, "mean" )

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

Si queremos probar el ajuste a un modelo unimodal (parábola), solo tenemos que añadir en la ecuación el término cuadrático de la variable cuantitativa (estadísticamente no significativo en este caso):

lm( cob ~ hum + I( hum ^ 2 ), data = ejemplo ) -> m3
summary( m3 )
plot( predictorEffects( m3 ) )

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

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

Los 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 ) -> m5
summary( m5 )
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 ) -> m6
summary( m6 )
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( m6, ~ hum ), multiline = TRUE, confint = list( style = "bands" ) )

Una función muy útil es anova (no debe confundirse con aov), que permite comparar dos o más modelos:

anova( m5, m6 )

Modelos lineales generalizados

Utilizaremos una matriz de datos simulados sobre la ocupación y reproducción de una rapaz territorial:

territorios
str( territorios )

En esta tabla aparecen símbolos NA (Not Availaible) en la variable pollos que representan valores perdidos (missing values), y que se corresponden con los casos en los que no hay ocupación. Con este ejemplo se explica bien la diferencia entre un 0 (la pareja no se reproduce o fracasa en la reproducción) y un NA (no hay valor de reproducción porque el territorio no está ocupado).

La ocupación, expresada con datos binarios (1/0), la analizaremos con modelos de regresión logística (distribución binomial de errores), usando la función glm:

glm( ocup ~ alt, family = binomial, data = territorios ) -> m7
summary( m7 )

Comprueba cómo se obtiene la curva ajustada utilizando la ecuación del modelo logístico y la función plogis, que realiza la transformación logística inversa:

\[\ln \left( \dfrac{ p }{ 1-p } \right)= b_0 + b_1 x \quad \rightarrow \quad p = \dfrac{ e^{ b_0 + b_1 x } }{ 1 - e^{ b_0 + b_1 x } }\]

plot( ocup ~ alt, data = territorios )
plogis( -9.180380 + 0.010784 * ( 700:1200 ) ) -> p
lines( 700:1200 , p )

La representación con los intervalos de confianza de este último modelo requiere un nuevo argumento (type = "response"):

plot( predictorEffects( m7 ), type = "response", ylim = c( 0, 1 ) )

Por su parte, la reproducción (productividad, expresada como número de pollos) la analizaremos con modelos de regresión discreta o de Poisson (datos de conteos con distribución de errores de Poisson):

glm( pollos ~ alt, family = poisson, data = territorios ) -> m8
summary( m8 )

Comprueba cómo se obtiene la curva ajustada utilizando la ecuación del modelo de Poisson:

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

plot( pollos ~ alt, data = territorios )
exp( 0.8271183 -0.0007943 * ( 700:1200 ) ) -> Ey
lines (700:1200, Ey )

La representación con los intervalos de confianza la obtendríamos con:

plot( predictorEffects( m8 ), type = "response" )

Interpretemos ahora los modelos con variables cualitativas:

glm( pollos ~ orient, family = poisson, data = territorios ) -> m9
summary( m9 )
plot( predictorEffects( m9 ), type = "response" )

Las diferencias entre los grupos de la variable cualitativa pueden analizarse modificando la modalidad de referencia con la función relevel. Por ejemplo:

glm( pollos ~ relevel( orient, 2 ), family = poisson, data = territorios ) -> m10
summary( m10 )

Utilizaremos también la función anova para testar la significación estadística del factor (variable orient):

anova( m10, test = "Chisq" )

En estos modelos con variables independientes cualitativas puede resultar interesante aplicar la opción de eliminar la estimación de la constante (~ -1), lo que facilita la interpretación de los coeficientes de cada modalidad (aunque no puedan evaluarse las diferencias estadísticas entre ellas):

glm( pollos ~ -1 + orient, family = poisson, data = territorios ) -> m11
summary( m11 )

Finalmente, el análisis “completo” de la ocupación y la reproducción de la rapaz lo obtendríamos introduciendo los términos de interacción entre las dos variables ambientales:

glm( ocup ~ alt * orient, family = binomial, data = territorios ) -> m12
anova( m12, test = "Chisq" )
plot( predictorEffects( m12, ~ alt ), multiline = TRUE, type = "response", confint = list( style = "bands" ) )
x11()
glm(pollos ~ alt * orient, family = poisson, data = territorios ) -> m13
anova( m13, test = "Chisq" )
plot( predictorEffects( m13, ~ alt ), multiline = TRUE, type= "response", confint = list( style = "bands" ) )

Selección de modelos e inferencia multimodelo

Estos procedimientos constituyen un nuevo paradigma para el análisis estadístico. En síntesis, un procedimiento de selección de modelos evalúa el ajuste a los datos de un conjunto de modelos “candidatos”, y realiza la estimación de parámetros en función de las estimaciones ponderadas de dicho parámetro en los diferentes modelos. Esta evaluación se realiza generalmente mediante el criterio de información de Akaike (AIC), o su versión corregida (AICc). Se busca la “mejor” explicación posible para los datos observados (best approximating model) y, bajo esta perspectiva, el uso de los valores de P de la estadística frecuentista no tiene sentido.

Usaremos de nuevo los datos del objeto territorios, para los que construiremos cinco modelos candidatos distintos, incluido el modelo nulo.

glm( ocup ~ 1, family = binomial, data = territorios ) -> mc0
glm( ocup ~ alt, family = binomial, data = territorios ) -> mc1
glm( ocup ~ orient, family = binomial, data = territorios ) -> mc2
glm( ocup ~ alt + orient, family = binomial, data = territorios ) -> mc3
glm( ocup ~ alt * orient, family = binomial, data = territorios ) -> mc4

Con el conjunto de modelos “candidatos”, usando la función model.sel del paquete MuMIn, obtendremos la tabla de selección:

model.sel( mc0, mc1, mc2, mc3, mc4 ) -> ms
ms
coef( ms )

Con la función de inferencia multimodelo (model.avg) estimaremos los diferentes coeficientes:

model.avg( ms ) -> ma
summary( ma )

La función proporciona dos tablas de coeficientes, una en la que se consideran todos los modelos candidatos (full average), y otra en la que se consideran solo los modelos donde aparece la variable en cuestión (conditional average). También resulta muy interesante la información sobre la “importancia relativa” de las variables (suma de los pesos de Akaike de los modelos en los que aparece cada variable).

sw( ma )

Los valores estimados por el modelo promediado, se obtienen con la función predict:

predict( ma, se.fit = TRUE, type = "response" ) -> map
map

La representación gráfica no se puede realizar con el paquete effects.

Evaluación

Realiza la tarea de evaluación de la práctica publicada en el 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( "territorios" ).