load( url( "https://webs.um.es/jfcalvo/mbc.RData" ) )
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
Actualizado el 14/11/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 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.
2. Desarrollo de la sesión
(A) 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.
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:
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 datos y funciones cargadas. 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 el programa.
Para realizar esta práctica se necesita una librería de R instalada por defecto (nlme), y otras que requieren una instalación adicional (effects, lme4 y MuMIn). Si no los tenemos instalados usaremos:
install.packages( c( "effects", "lme4", "MuMIn" ) )Una vez instalados los cargaremos todos en memoria:
library( nlme )
library( effects )
library( lme4 )
library( MuMIn )(B) ANOVA, regresión lineal y ANCOVA
Comenzaremos realizando análisis de la varianza (ANOVA) utilizando los datos simulados contenidos en el objeto ejemplo:
ejemploPodemos 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( ejemplo$cov, ejemplo$y, col = ejemplo$trat )O mejor:
plot( y ~ cov, data = ejemplo, col = trat, 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( y ~ trat, data = ejemplo ) -> m1
summary( m1 )
model.tables( m1, "mean" )El análisis anterior es equivalente a una prueba de la \(t\) de Student:
t.test( y ~ trat, var.equal = TRUE, data = ejemplo )Cuando la variable independiente es cuantitativa haremos una regresión (modelo lineal, con errores de distribución normal o gaussiana):
lm( y ~ cov, data = ejemplo ) -> m2
summary( m2 )
plot( y ~ cov, data = ejemplo )
abline( 14.5833, 2.9985 )Una representación adecuada del modelo debe incluir el intervalo de confianza. Para ello utilizaremos la función predictorEffects del paquete effects:
plot( predictorEffects( m2 ) )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( y ~ cov + I( cov ^ 2), data = ejemplo ) -> m3
summary( m3 )
plot( predictorEffects( m3 ) )También podemos usar un modelo lineal para hacer un análisis equivalente al ANOVA:
lm(y ~ trat, 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):
lm( y ~ trat + cov, data = ejemplo ) -> m5
summary( m5 )
plot( y ~ cov, col = trat, data = ejemplo )
abline( 21.4583, 2.9985 )
abline( 21.4583 -13.75, 2.9985, col = 2 )Interpretemos un modelo más complejo que incluya términos de interacción entre las variables:
lm(y ~ trat * cov, data = ejemplo ) -> m6
summary( m6 )
abline( 31.3667, 1.197, lwd = 2)
abline( 31.3667 -33.5667, 1.197 + 3.603, lwd = 2, col = 2 )La representación gráfica de este último modelo, con los intervalos de confianza, la haríamos con:
plot( predictorEffects( m6, ~ cov ), 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 )(C) 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 Available) 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( -8.869028 + 0.010527 * ( 700:1100 ) ) -> p
lines( 700:1100 , 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.8613412 -0.0008336 * ( 700:1100 ) ) -> Ey
lines (700:1100, Ey )
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" )Utilizaremos también la función anova para testar cada variable/factor:
anova( m9, test = "Chisq" )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 )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" ) )(D) Modelos mixtos
Los modelos de regresión mixtos permiten analizar de manera apropiada la falta de independencia en datos que se agrupan en función de individuos o unidades experimentales o de muestreo y que corresponden habitualmente a medidas u observaciones repetidas. En un modelo mixto se asume que elemento de agrupación, denominado factor (o efecto) aleatorio es una muestra aleatoria de una población mayor. Por otra parte, las variables potencialmente explicativas de la respuesta biológica objeto de la investigación se denominan factores (o efectos) fijos.
Retomaremos los datos del objeto ejemplo, analizados al inicio de la práctica obviando la falta de independencia, y que serán analizados aquí con modelos más apropiados a su estructura. En esta matriz variable ind representa el factor aleatorio y las variables trat y cov se consideran efectos fijos.
ejemploVisualicemos inicialmente la estructura agrupada de los datos en una gráfica:
plot( groupedData( y ~ cov | ind, data = ejemplo ) )Usaremos la función lme para ejecutar un modelo mixto sencillo:
lme( y ~ cov, random = ~ 1 | ind, data = ejemplo ) -> m14
summary( m14 )El análisis proporciona una única respuesta al efecto fijo para los cuatro individuos, y otra distinta para cada uno de ellos (efecto aleatorio), aunque en este caso todas con la misma pendiente:
m14$fit
plot( m14$fit[ , 2 ] ~ cov, data = ejemplo, col = ind )
points( m14$fit[ , 1 ] ~ cov, data = ejemplo, pch = 19 )Un modelo más complejo podría contemplar la correlación entre constante y pendiente:
lme( y ~ cov, random = ~ 1 + cov | ind, data = ejemplo ) -> m15
summary( m15 )
plot( m15$fit[ , 2 ] ~ cov, data = ejemplo, col = ind )
points( m15$fit[ , 1 ] ~ cov, data = ejemplo, pch = 19 )Con la función anova podemos comprobar que el segundo modelo mejora considerablemente el ajuste:
anova( m14, m15 )Para ejecutar modelos lineales mixtos generalizados usaríamos la función glmer del paquete lme4. Como ejemplo analizaremos los datos de ocupación del objeto territorios, considerando la variable terr como factor aleatorio (1 | terr) y la variable ambiental alt como factor fijo. Observa que con glmer la forma de especificar los factores aleatorios es diferente a la de lme:
glmer( ocup ~ alt + ( 1 | terr ), family = binomial, data = territorios ) -> m16
summary( m16 )En este caso la varianza entre territorios es 0, por lo que el modelo mixto apenas difiere del “clásico” (modelo m7):
summary ( m7 )
anova( m16, m7 )(E) 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 ) -> mc4Con 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 )La representación gráfica no se puede realizar con el paquete effects, por lo que hay que recurrir a los valores estimados por el modelo promediado, que se obtienen con la función predict:
predict( ma, se.fit = TRUE ) -> map
map
plot( ocup ~ alt, data = territorios, pch = 19 )
lines( plogis( map$fit )[ order( alt ) ] ~ sort( alt ), data = territorios, col = 4, lwd = 2 )
lines( plogis( map$fit - map$se * 1.96 )[ order( alt ) ] ~ sort( alt ), data = territorios, col = 4 )
lines( plogis( map$fit + map$se * 1.96 )[ order( alt ) ] ~ sort( alt ), data = territorios, col = 4 )3. Evaluación
Realiza la tarea de evaluación de la práctica publicada en el 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
Utiliza la función info. Por ejemplo: info("territorios").