Paquetes R

library(drc)
library(primer)
library(readxl)
library(stats)
library(kableExtra)

Modelos Continuos

Cuando los eventos estocásticos no son los principales determinantes de los parámetros demográficos, generalmente podemos describir el crecimiento (o disminución) poblacional, utilizando modelos empíricos, basados en funciones matemáticas continuas, usualmente una función del tiempo.

Si logramos obtener los parámetros de una función que se ajuste a los datos históricos de una población, podemos utilizar la misma para predecir las variaciones poblacionales en función del tiempo, y simular el comportamiento de la población al modificar los parámetros, y asi estudiar el efecto de nuevos escenarios e inclusive introducir efectos aleatorios.

Conocemos modelos clásicos (exponencial y logístico) que ayudan a entender el crecimiento poblacional de forma general, y que en algunos casos se ajustan a poblaciones, generalmente en condiciones controladas o estables (crecimiento bacteriano); otra estrategia para encontrar modelos que se ajusten a los datos, consiste en usar métodos estadísticos, como la regresión lineal, pero este procedimiento ajusta a una línea recta o datos que se puedan linearizar, mediante transformaciones.

Ejemplo 1

En este ejercicio intentaremos ajustar datos poblacionales a diversas funciones matemáticas (diferentes a las clásicas). El ejercicio está basado en un artículo de Szabelska y col. (2010), sobre comparación de modelos de crecimiento usando R. La meta es encontrar la función (modelo) que describe mejor sus datos (bajo AIC), con coeficientes significativamente diferentes de 0. A partir de este modelo podemos simular y proyectar el crecimiento poblacional.

#data.frame de prueba
dataset1 <- read_excel("dataset1.xlsx")
#gráfica de los puntos
plot(dataset1$Time1, dataset1$Response1, xlab = "Tiempo", ylab = "N")

Figura 1. Gráfica de los datos para ser ajustados a un modelo de crecimiento.

Usaremos la función drm del paquete drc. El parámetro fct sirve para indicar el modelo a usar.

tablamod <- data.frame(
  Modelo = c("Exponencial", "Gompertz", "Logístico", "Log-logístico", "Weibull"),
  Parametro.fct = c("EXD.3()", "G.4()", "L.5()", "LL.5()", "W1.4()")
)
kable(tablamod, caption = "Tabla 1.  Códigos para modelos en fct") %>%
  kable_styling(full_width = F) %>%
  column_spec(1, bold = T, border_right = T) %>%
  column_spec(2, width = "10em")
Tabla 1. Códigos para modelos en fct
Modelo Parametro.fct
Exponencial EXD.3()
Gompertz G.4()
Logístico L.5()
Log-logístico LL.5()
Weibull W1.4() 
#cálculo de parámetros de cada modelo (repetir c/u) - aqui con Gompertz
result.G <- drm(dataset1$Response1~dataset1$Time1, data = dataset1, fct = G.4())
#parámetros y prueba estadística
summary(result.G)
## 
## Model fitted: Gompertz (4 parms)
## 
## Parameter estimates:
## 
##                  Estimate  Std. Error  t-value   p-value    
## b:(Intercept)  -0.0209818   0.0010291 -20.3892 < 2.2e-16 ***
## c:(Intercept)   8.9892304   2.0628820   4.3576 0.0001831 ***
## d:(Intercept) 240.7688604   3.2933892  73.1067 < 2.2e-16 ***
## e:(Intercept) 135.2307582   1.4890631  90.8160 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  4.891884 (26 degrees of freedom)
#gráfica
plot(result.G, xlab = "Tiempo", ylab = "N")

Figura 2. Gráfica de los datos y el modelo calculado con la función drm.

Cuando se usan varios modelos, y todos presentan parámetros significativos, se puede usar el estadístico AIC para seleccionar el mejor. Un valor bajo de este estadístico, significa un mejor ajuste del modelo a los puntos.

#evaluación con AIC
aic <- AIC(result.G)
sprintf("AIC: %.1f", aic)
## [1] "AIC: 186.1"

Modelos con Matrices

El ciclo de vida de un organismo se puede representar gráficamente utilizando un diagrama de clases de edades o etapas y los ‘flujos’ o transiciones entre ellos. tomado de Geib & Galen (2012)

Figura 3. Diagrama del ciclo de vida de dos especies de Trifolium.

Para realizar el análisis de un ciclo de vida que cuente con los parámetros demográficos, podemos utilizar un modelo basado en una matriz de transición. Una vez que se obtenga la matriz, podemos proyectar el crecimiento poblacional (no denso-dependiente), y simular escenarios con variaciones de las transiciones entre etapas. Los datos para construir el modelo (transiciones, fertilidad), y las proyecciones son discretos, basados en iteraciones temporales.

Ejemplo 2

Vamos a proyectar la población usando datos de la Palma de Sierra (Prestoea acuminata), tomados en la Estación Biológica El Verde. Estaremos asumiendo que las transiciones y el parámetro \(\lambda\) se calcularon para una población de estructura estable.

Vamos a hacer proyecciones poblacionales cambiando las condiciones iniciales de la población, y simulando el escenario luego del paso del huracán María.

Construcción de matriz:

A <- matrix(c(0, 0, 0, 28, 
              0.1, 0.3, 0, 0, 
              0, 0.3, 0.1, 0, 
              0, 0, 0.2, 0.8), 
            nr=4, byrow=TRUE)
sprintf("Matriz de Transición")
## [1] "Matriz de Transición"
A
##      [,1] [,2] [,3] [,4]
## [1,]  0.0  0.0  0.0 28.0
## [2,]  0.1  0.3  0.0  0.0
## [3,]  0.0  0.3  0.1  0.0
## [4,]  0.0  0.0  0.2  0.8

Ahora creamos el vector de población inicial (\(N_0\)), y proyectamos la población multiplicando por la matriz de transición:

N0 <- matrix(c(875, 50, 50, 25), ncol=1)
sprintf("Vector población inicial N0")
## [1] "Vector población inicial N0"
N0
##      [,1]
## [1,]  875
## [2,]   50
## [3,]   50
## [4,]   25
N1 <- A %*% N0
sprintf("Vector población a t=1")
## [1] "Vector población a t=1"
N1
##       [,1]
## [1,] 700.0
## [2,] 102.5
## [3,]  20.0
## [4,]  30.0

Podemos hacer la proyección a varios años y graficarla:

years <- 30
N.projections <- matrix(0, nrow=nrow(A), ncol=years+1) 
N.projections[,1] <- N0
  for(i in 1:years) N.projections[,i+1] <- A%*%N.projections[,i]
matplot(0:years, t(N.projections), type="l", lty=1, col=1:4,
        ylab="Stage Abundance", xlab="Year")
legend('topleft', legend=c("S1","S2", "S3", "S4"), lty=1, 
       col=1:4, bty='n')

Figura 4. Proyección poblacional para las etapas (‘stages’) del ciclo de vida de la población.

A continuación calculamos \(\lambda\) para cada iteración (anual, usualmente) y la graficamos:

N.totals <- apply(N.projections, 2, sum)
Rs <- N.totals[-1] / N.totals[-(years+1)]
plot(0:(years-1), Rs, type="b", xlab="año", ylab="lambda")

Figura 5. Gráfica del parámetro \(\lambda\) para cada iteración.

Hemos realizado proyecciones utilizando el modelo de matriz de transición, y observamos que \(\lambda\) alcanzó un valor estable. Si este es el caso, ocurrirá cuando la población tenga una estructura de edades o etapas estable (proporciones constantes). Podemos obtener estos valores finales de \(\lambda\) y el vector de estructura estable, calculando los autovalores (‘eigenvalues’).

eigs.A <- eigen(A)
dom.pos <- which.max( eigs.A[["values"]] ) 
Lambda <- Re(eigs.A[["values"]][dom.pos])
sprintf("Lambda")
## [1] "Lambda"
Lambda
## [1] 1.035663
w <- Re(eigs.A[["vectors"]][,dom.pos])
ssd <- w/sum(w)
sprintf("Proporción de Etapas Estable")
## [1] "Proporción de Etapas Estable"
round(ssd, 3)
## [1] 0.822 0.112 0.036 0.030

Otra cantidad importante en el análisis del modelo de matriz de transición es el cálculo del valor reproductivo, una medida de la capacidad de cada etapa de aportar a la reproducción en la población.

M <- eigen(t(A))
v <- Re(M$vectors[,which.max(Re(M$values))])
RV <- v/v[1]
sprintf("Valores Reproductivos")
## [1] "Valores Reproductivos"
RV
## [1]   1.00000  10.35663  25.39665 118.81357

Para responder a la pregunta: ¿cuál transición aporta más al valor de \(\lambda\)?, se realiza un análisis de sensibilidad.

vw.s <- v %*% t(w)
S <- round(vw.s/as.numeric(v%*%w), 3)
sprintf("Matriz de Sensibilidad")
## [1] "Matriz de Sensibilidad"
S
##        [,1]  [,2]  [,3]  [,4]
## [1,]  0.126 0.017 0.006 0.005
## [2,]  1.309 0.178 0.057 0.048
## [3,]  3.211 0.436 0.140 0.119
## [4,] 15.022 2.042 0.655 0.556

Igualmente podemos responder a la pregunta: ¿cuál cambio en una transición afecta más a \(\lambda\)?; para esto calculamos la matriz de elasticidad.

elas <- (A/Lambda) * S
sprintf("Matriz de Elasticidad")
## [1] "Matriz de Elasticidad"
round(elas, 3)
##       [,1]  [,2]  [,3]  [,4]
## [1,] 0.000 0.000 0.000 0.135
## [2,] 0.126 0.052 0.000 0.000
## [3,] 0.000 0.126 0.014 0.000
## [4,] 0.000 0.000 0.126 0.429

Referencias

Geib, J.C., Galen, C. 2012. Tracing impacts of partner abundance in facultative pollination mutualisms: from individuals to populations. Ecology 93:1581–1592. https://doi.org/10.1890/11-1271.1

Stevens, M. H. H. 2009. A Primer of Ecology with R, First Edition. ed. Springer, US.

Szabelska A., Siatkowski M., Goszczurna T., Zyprych J. 2010. Comparison of growth models in package R. Nauka. Przyr. Technol. 4(4):1-9, #50.

EJERCICIOS

  1. Abundancia de salmones en río
library(readxl)
salmon <- read_excel("salmon_abundance.xlsx", 
    sheet = "abundance")
#gráfica de los puntos
plot(salmon$year, salmon$abund, xlab = "Tiempo", ylab = "N")

library(drc)
#cálculo de parámetros de cada modelo (repetir c/u)
result.M <- drm(salmon$abund~salmon$year, data = salmon, fct = G.4())
#parámetros y prueba estadística
summary(result.M)
## 
## Model fitted: Gompertz (4 parms)
## 
## Parameter estimates:
## 
##                  Estimate  Std. Error     t-value p-value    
## b:(Intercept) -4.5324e-01  1.0000e+01     -0.0453  0.9645    
## c:(Intercept)  3.5503e+04  9.9999e+00   3550.2933  <2e-16 ***
## d:(Intercept) -3.8675e+05  1.0000e+01 -38675.2028  <2e-16 ***
## e:(Intercept)  2.0760e+03  1.0000e+01    207.6036  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  12263.14 (14 degrees of freedom)
#gráfica
plot(result.M, xlab = "Tiempo", ylab = "N")

aic <- AIC(result.M)
sprintf("AIC: %.1f", aic)
## [1] "AIC: 395.5"
  1. Proyección poblacional de orquídea