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 aplicación de varias técnicas de análisis de regresión. Utilizaremos métodos 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 online Análisis estadístico de datos biológicos en R, donde encontrarás más información, ejercicios, ejemplos y métodos de análisis más avanzados.

Preparación

Para la realización de la 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 online Primeros pasos en R.

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. Puedes utilizar 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://webs.um.es/jfcalvo/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 R.

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

El gráfico anterior puede mejorarse con la introducción de diversos argumentos. Por ejemplo:

plot( predictorEffects( m6, ~ hum ), multiline = TRUE, confint = list( style = "bands" ),
      main = "", xlab = "Humedad del suelo", ylab = "Cobertura (%)",
      lattice = list( key.args = list( x = 0.7, y = 0.2, title ="Tipo de hábitat" ) ) )
points( cob ~ hum, data = ejemplo, col = hab, pch = 19, cex = 1.5 )

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

Para la realización de ejercicios con modelos lineales generalizados utilizaremos una matriz de datos simulados sobre la ocupación y reproducción de una rapaz territorial en 70 territorios localizados a diferentes altitudes y orientaciones:

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 la variable cualitativa:

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( m9, 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. Para la ocupación:

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

Y para la reproducción:

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

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 el seminario online Análisis estadístico de datos biológicos.

Selección de modelos e inferencia multimodelo

Los procedimientos de selección de modelos y la inferencia multimodelo 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, que se calculan como media ponderada usando los pesos de Akaike (\(w\)):

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). Por lo general, se recomienda utilizar los coeficientes de la tabla full 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

Aunque la representación gráfica no puede realizarse con el paquete effects, podemos utilizar la función plot.ic (disponible en el archivo mbc.RData), aplicada sobre una tabla con los valores originales de las variables, los valores estimados por predict (map$fit, en este caso), y los valores de los límites inferior y superior del intervalo de confianza (calculados con map$se.fit).

data.frame(territorios, map$fit, map$fit - 1.96 * map$se.fit, map$fit + 1.96 * map$se.fit) -> map2
map2

Como el modelo predice una curva para cada orientación, tendremos que seleccionar los casos (filas). Por ejemplo, para la orientación Norte:

subset( map2, map2$orient == "Norte" ) -> norte
plot.ic( x = norte$alt, y = norte[ , 5 : 7 ], xlab = "Altitud", ylab = "Probabilidad de ocupación" )

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