load( url( "http://webs.um.es/jfcalvo/mbc.RData" ) )
Análisis de viabilidad de poblaciones
Práctica 5, Métodos en Biología de la Conservación
Máster en Áreas Protegidas, Recursos Naturales y Biodiversidad
Actualizado el 20/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 analizaremos las características de los modelos estocásticos de crecimiento densoindependiente para realizar análisis de viabilidad de poblaciones basados en censos (series temporales), aprendiendo a interpretar el significado de la varianza de las tasas de crecimiento y su influencia en las estimas de la probabilidad o riesgo de cuasiextinción. Realizaremos pruebas de diagnóstico sobre las asunciones de los métodos aplicados con modelos densoindependientes (estabilidad temporal de los parámetros, ausencia de autocorrelación y de outliers), para centrarnos posteriormente en las características de los modelos estocásticos de crecimiento densodependiente y los efectos de la autocorrelación positiva de las tasas de crecimiento poblacionales en la probabilidad de extinción.
En la segunda parte de la práctica repasaremos inicialmente las características de las matrices de proyección y de los análisis demográficos de las mismas, con objeto de realizar análisis de viabilidad de poblaciones demográfico. Finalmente realizaremos análisis de sensibilidad y elasticidad, evaluando el efecto en la tasa de crecimiento anual de cambios en los elementos de las matrices de proyección.
2. Desarrollo de la sesión
(A) Preparación
Una vez iniciado R cargaremos el archivo de datos y funciones (mbc.RData) desde el servidor:
Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "condor" ) o info( "extCDF" ). Teclea info() para más información.
Para el desarrollo de uno de los ejercicios incluidos en este guion también se necesita el paquete MuMIn. Si no lo tenemos instalado utilizaremos:
install.packages( "MuMIn" )Una vez instalado lo cargaremos en memoria:
library( MuMIn )(B) Modelos estocásticos
Los modelos estocásticos incorporan variabilidad en los parámetros demográficos, por lo que, a diferencia de los modelos deterministas, la estimación del tamaño poblacional en el futuro se basa en criterios probabilísticos. Una aplicación de gran interés de estos modelos es la estimación de probabilidades de extinción, o más concretamente de cuasiextinción (la probabilidad de alcanzar un tamaño de población crítico \(N_{x}\)). La función extCDF puede simular modelos estocásticos, lo cual nos permitirá comprobar el efecto de la varianza en las probabilidades de cuasiextinción. Utilizaremos el argumento ejemplo = TRUE, para realizar una primera simulación:
extCDF( ejemplo = TRUE )Sin utilizar por el momento los parámetros K, theta y rho, comprobaremos que la probabilidad de extinción aumenta conforme aumenta el valor de la varianza, por ejemplo:
extCDF( ejemplo = TRUE, s2 = 0.15 )Dado que el número de simulaciones por defecto es bajo (20), las probabilidades estimadas son muy variables. Necesitaríamos un número de simulaciones mucho mayor para obtener una estima más precisa (por ejemplo 10000):
extCDF( ejemplo = TRUE, s2 = 0.15, sim = 10000 )(C) Ejemplo de estimación de parámetros de un modelo estocástico
Utilizaremos un sencillo ejemplo de una serie temporal de cuatro años. En un modelo estocástico, \(\lambda\) es variable, de manera que:
\[N_{4} = \lambda_{1}\ \lambda_{2}\ \lambda_{3}\ \lambda_{4}\ N_{0}\]
Para esta serie, consideraremos un tamaño inicial N = 100, y cuatro valores distintos de \(\lambda_t\) (0.7, 1.1, 0.9 y 1.2):
N0 <- 100
lambdas <- c( 0.7, 1.1, 0.9, 1.2 )0.7 * N0 -> N1
1.1 * N1 -> N2
0.9 * N2 -> N3
1.2 * N3 -> N4
N0; N1; N2; N3; N4Lo que es equivalente a:
0.7 * 1.1 * 0.9 * 1.2 * N0o también:
prod( lambdas ) * N0Es interesante comprobar que la media aritmética de las \(\lambda_t\) no describe adecuadamente la dinámica de la población.
mean( lambdas )
mean( lambdas ) ^ 4 * N0Lo correcto es utilizar la media geométrica:
\[\lambda_{G} = \left( \prod_{}^{}\lambda_{t} \right)^{\left( \frac{1}{t} \right)}\]
que en R calculamos así:
prod( lambdas ) ^ ( 1 / 4 )Ahora sí podemos comprobar que con esta media geométrica obtenemos el tamaño de población esperado en 4 años:
\[N_{t} = \lambda_{G}^{t}\ N_{0}\]
0.9549456 ^ 4 * N0Comprueba que este valor de \(\lambda_G\) se obtiene también, de forma más simple, calculando:
(N4 / N0) ^ ( 1 / 4 )Cuando se trabaja con modelos estocásticos se utilizan habitualmente los logaritmos de las \(\lambda_t\), cuya media (aritmética en este caso), es el parámetro \(\mu\) (equivalente en un modelo estocástico a la \(r\) de un modelo determinista). Comprueba que \(\mu\) puede calcularse también como el logaritmo de la media geométrica de las \(\lambda_t\):
\[\mu = \frac{\sum_{}^{}{\log\lambda_{t}}}{t} = \log\lambda_{G}\]
En R:
mean( log( lambdas ) )
log( prod( lambdas ) ^ ( 1 / 4 ) )El parámetro de varianza del modelo estocástico (\(s^2\)) corresponde precisamente a la varianza de \(\mu\), que se calcula como:
var( log( lambdas ) )(D) Ejemplo del cóndor californiano
Utilizaremos los datos de una serie temporal correspondiente al cóndor californiano (Gymnogyps californianus):
condorLa tabla contiene dos variables: el año (year) y el tamaño de la población (Nt). La información adicional sobre los datos cargados puede consultarse usando: info( "condor " ).
Empezaremos a analizar la serie temporal del cóndor, para obtener los parámetros \(\mu_t\) y \(s^2\). En primer lugar calcularemos las \(\lambda_t\) para cada año:
lambdas <- condor$Nt[ 2:16 ] / condor$Nt[ 1:15 ]
lambdasComprueba, por ejemplo, que \(\lambda_6\) (es decir, \(\lambda_{1970}\)) es el cociente entre el tamaño poblacional de 1971 y el de 1970.
A continuación, a partir del vector de lambdas, calcularemos su media geométrica \(\lambda_G\) y finalmente \(\mu\) y \(s^2\).
prod( lambdas ) ^ ( 1 / 15 )
log( prod( lambdas ) ^ ( 1 / 15 ) )
mean(log( lambdas ) )
var( log( lambdas ) )Cuando la serie temporal de censos no está completa, este procedimiento de cálculo no puede aplicarse, y hay que utilizar otro método algo más complejo (detalles en Morris & Doak 2002). La función extCDF permite realizar estos cálculos de manera automática, con sus intervalos de confianza al 95%:
extCDF( condor )Cuando se proporciona una tabla de datos (con dos columnas: year y Nt), la función extCDF también proporciona la función de distribución acumulada (CDF), para un valor por defecto de \(N_{x} =\) 1 (que puede modificarse), considerando un tamaño inicial de población \(N_{c}\) igual al censo del último año de la tabla. En este caso, extCDF realiza una estimación analítica de la CDF con intervalos de confianza bootstrap (nboot = 500, level = 0.95), basada en la aproximación de difusión (Morris & Doak 2002).
Si quisiéramos realizar un análisis mediante simulación, introduciríamos los parámetros directamente en la función. Podemos abrir una nueva ventana gráfica para comparar ambas curvas:
x11()
extCDF( mu = -0.0768453, s2 = 0.1175639, Nc = 12, Nx = 1)La función extCDF realiza por defecto tres repeticiones (tres estimaciones de la CDF, rep = 3) de 50000 trayectorias poblacionales aleatorias simuladas (sim = 50000) cada una, para un horizonte temporal de 100 años (tiempo = 100).
(E) Diagnósticos
La validez de los análisis de viabilidad de poblaciones basados en censos depende de que se cumplan una serie de asunciones relacionadas con: (1) la estabilidad temporal de los parámetros poblacionales, (2) la ausencia de autocorrelación interanual entre dichos parámetros y (3) la ausencia de outliers, es decir, datos que se “salen de lo normal” y que ejercen una influencia considerable sobre las estimas de viabilidad. Para este apartado utilizaremos la serie temporal del oso grizzly (Ursus arctos horribilis) de Yellowstone:
grizzlyUtilizaremos la función extCDF con el argumento tests = TRUE. Esta opción proporciona pruebas estadísticas de correlación entre los valores \(\log\lambda_{t}\) y los \(N_{x}\) (densodependencia), autocorrelación (\(\rho\)) de los residuos y diversos diagnósticos de regresión (medidas de influencia).
extCDF( grizzly, tests = TRUE )Las medidas de influencia se estiman con la función influence.measures, que proporciona diversos estadísticos que miden la influencia de una observación en el modelo: DFBETAS, DFFITS, covariance ratios, distancias de Cook y los elementos de la diagonal de la matriz \(\mathbf{H}\) o hat matrix (también conocidos como leverages). Las observaciones o casos (años) que son influyentes según cualquiera de dichas medidas se marcan con un asterisco. En el caso del grizzly no se aprecia densodependencia ni autocorrelación, pero sí la presencia de un outlier (observación número 25), que corresponde al valor excepcionalmente alto de \(\lambda\) entre los años 1983 y 1984. Podría valorarse, por tanto, la posibilidad de eliminar ese valor en los datos para realizar los análisis de viabilidad.
(F) Modelos densodependientes
Analizaremos el caso de una población en la que se observa una relación negativa estadísticamente significativa entre \(\log\lambda_{t}\) y \(N_{t}\), lo que requiere el ajuste de los datos a un modelo densodependiente. Probaremos con el modelo de Ricker y el modelo \(\theta\)-logístico.
\[N_{t + 1} = N_{t}e^{r\left( 1 - \frac{N_{t}}{K} \right) + \varepsilon}\]
\[N_{t + 1} = N_{t}e^{r\left( 1 - \left( \frac{N_{t}}{K} \right)^{\theta} \right) + \varepsilon}\]
Utilizaremos los datos temporales de la población “JRC” de la mariposa Euphydryas editha bayensis en la Bahía de San Francisco (1960-1986). La revisión de los diagnósticos muestra la existencia de densodependencia negativa:
extCDF( JRC, tests = TRUE )Ajustaremos los datos a tres modelos (densoindependiente, Ricker y \(\theta\)-logístico), estimaremos los parámetros de cada modelo (\(r\), \(K\), \(\theta\)) y la varianza residual (\(s^{2}\)) y los compararemos mediante el criterio de información de Akaike corregido (AICc):
logL <- extCDF( JRC )$logL
N <- JRC$Nt[ -length( JRC$Nt ) ]
cbind( logL, N)m1 <- lm( logL ~ 1 )
coef( m1 )
var( summary( m1 )$res )m2 <- nls( logL ~ r * ( 1 - ( N / K ) ) )
coef( m2 )
var( summary( m2 )$res )m3 <- nls( logL ~ r * ( 1 - ( N / K ) ^ theta ) )
coef(m3)
var( summary( m3 )$res )AICc( m1 ); AICc( m2 ); AICc( m3 )Para estimación de la probabilidad de extinción de la población utilizaremos de nuevo la función extCDF, considerando los parámetros estimados por el mejor modelo y un umbral de cuasiextinción de 20 individuos:
extCDF( r = 0.3458, s2 = 1.1151, K = 846, Nc = 94, Nx = 20 )Además de los parámetros \(r\), \(s^{2}\) y \(K\), la función extCDF admite valores de los parámetros \(\theta\) (theta) y \(\rho\) (rho), para aquellos casos en los que el mejor modelo sea el theta-logístico y en los que exista autocorrelación.
(G) Análisis de viabilidad demográficos
Los análisis demográficos de matrices de proyección los podemos realizar con la función matriz. Recordaremos las características de este tipo de análisis con la matriz de proyección de la ballena franca del Norte (Eubalaena glacialis):
matriz( ballena, plots = TRUE )Para realizar un AVP demográfico utilizaremos la función extMAT con cinco matrices de proyección anuales (2001/2002 a 2005/2006) de la población de dragoncillo de roca (Antirrhinum subbaeticum) en la zona Benizar-Hondares (Iriondo et al. 2009). Las matrices se encuentran en un objeto de tipo “lista” (asBH) en el archivo mbc.RData, y corresponden a una población estructurada en tres clases de individuos: plántula, vegetativo y reproductor. [También están disponibles las matrices de otra población de esta especie (Bogarra-Potiche-Mundo): asBPM.]
asBHAprovecharemos el ejercicio para comprobar que la estructura inicial de la población condiciona los resultados del análisis de viabilidad. Por ejemplo, considerando un umbral de cuasiextinción \(N_{x} =\) 20, y partiendo de un tamaño inicial de población de 45 individuos, pero con diferente estructura poblacional:
extMAT( asBH, n0 = c( 15, 15, 15) , Nx = 20 )extMAT( asBH, n0 = c( 25, 15, 5 ), Nx = 20 )(H) Análisis de sensibilidad
Para realizar análisis de sensibilidad utilizaremos como ejemplo la matriz de proyección de una población de ganso emperador (Anser canagicus) con cuatro clases de edad (0 años, 1 año, 2 años y 3 o más años). El argumento SE = T en la función matriz proporciona las matrices de sensibilidad y elasticidad, y un gráfico representando la simulación del efecto de cambios proporcionales en los elementos de la matriz en la tasa anual de crecimiento.
ganso
matriz( ganso, SE = TRUE )3. Evaluación
Practica con otros ejemplos utilizando las siguientes series de datos disponibles (con info se obtiene información sobre todas ellas): cactus, crane, elk, grizzly, JRH, parrot, perdicera, real, redstart, robin, turtle, turtleB, wolf. Elabora un informe de una de las dos alternativas propuestas:
Elige una serie para la cual la correlación entre \(\log\lambda_{t}\) y \(N_{t}\) sea negativa y estadísticamente significativa: estima los parámetros de un modelo densodependiente y realiza un AVP. Elige un umbral de cuasiextinción (\(N_{x}\)) de acuerdo con las características de la población.
Elige una serie para la que se detecte autocorrelación (\(\rho\)) positiva y estadísticamente significativa. Elige un umbral de cuasiextinción (\(N_{x}\)) de acuerdo con las características de la población y compara los resultados de dos AVP, considerando y no considerando el parámetro de autocorrelación.
[Recuerda que la función extCDF puede utilizarse para simular modelos densoindependientes usando K = 1e+20.]
4. Bibliografía
Iriondo JM, Albert MJ, Giménez Benavides L, Domínguez Lozano F, Escudero A (Eds.) (2009) Poblaciones en Peligro: Viabilidad Demográfica de la Flora Vascular Amenazada de España. Dirección General de Medio Natural y Política Forestal (Ministerio de Medio Ambiente, y Medio Rural y Marino), Madrid.
Morris WF, Doak DF (2002) Quantitative Conservation Biology. Sinauer, Sunderland, MA, USA.
Descripción de los datos y funciones
Utiliza la función info. Por ejemplo: info( "condor" ) o info( "extCDF" ).