Crecimiento poblacional con estructura de edades o de estados. Aplicación en R (paquete primer)

• Crecimiento con estructura de edades o de estados

###### 3) CRECIMIENTO DENSO-INDEPENDIENTE DEMOGRÁFICO (ESTRUCTURA ETARIA O DE ESTADOS)
#EJ construimos la matriz de proyección poblacional para una especie hipotética de planta perenne con semillas que duran un año.
M = matrix(1:4, nr = 2, byrow = T); N = matrix(c(10, 20, 30, 40), nr = 2); M %*% N
#crecimiento estado-estructurado en un paso
A= matrix(c(0, 0.5, 20, 0.3, 0, 0, 0, 0.5, 0.9), nr = 3, byrow = TRUE); N0 = matrix(c(100, 250, 50), ncol = 1); N1 = A %*% N0  #creamos la proyección poblacional y un vector de clases de estado para las abundancias para el año cero
#crecimiento estado-estructurado en múltiples pasos
years= 6; 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:3, col = 1, ylab = "Stage Abundance", xlab = "Year");    legend("topleft", legend = c("Seeds", "Small Adult", "Large Adult"), lty = 1:3, col = 1, bty = "n")
#tasa de crecimiento anual
N.totals= apply(N.projections, 2, sum);    Rs =N.totals[-1]/N.totals[-(years + 1)]
plot(0:(years - 1), Rs, type = "b", xlab = "Year", ylab = "R")
#analizando la matriz de proyección: análisis propio
eigs.A= eigen(A) #el valor propio dominante o primero, es la tasa finita asintótica a largo plazo del incremento lambda, sus correspondientes vectores propios proveen la distribución de estados estables. El vector propio dominante izquierdo provee los valores reproductivos.
dom.pos= which.max(eigs.A[["values"]]);L1 = Re(eigs.A[["values"]][dom.pos]) #encontrando tasa de crecimiento finita (lambda). Nos quedamos con la parte real utilizando Re.
t=20; Nt= N0/sum(N0); R.t= numeric(t); for (i in 1:t){ Nt1=A %*% Nt; R=sum(Nt1)/sum(Nt); Nt= Nt1/sum(Nt1); R.t[i]=R } #Otra forma de encontrar a lambda1 es iterar el crecimiento poblacional un gran número de veces, cuando t crece la tasa de crecimiento anual Nt+1/Nt se aproxima a lambda1.
par(mar = c(5, 4, 3, 2)); plot(1:t, R.t, type = "b", main = quote("Convergence Toward " * lambda)); points(t, L1, pch = 19, cex = 1.5) #comparamos los resultados de estimar lambda1.
#distribuciones estables de estados (abundancia relativa de los diferentse estados de historia de vida): sea A*w=lambda*w  entonces SSD=w1/sum(w1). Propiedad: si todas las tasas demográficas (elementos de la matriz de proyección) permanecen constantes, sus estructuras de estado se aproximará a una distribución de estados estable, donde el número relativo de individuos en casa estado es constante. Si la població no está creciendo (lambda=1) los parámetros demográficos permanecen constantes, entonces la población tiene distribuciónd e estados estacionaria, donde la abundancia absoluta y no la relativa cambia.
w=Re(eige.A[["vectors"]][,dom.pos]); ssd=w/sum(w);round(ssd,3) #SSD=w1/sum(w1), donde S es el número de estados y w1 el vector propio de lambda1. Caundo la población alcanza su distribución de estados estable, crece exponencialmente: Nt=A^t*No (representación matricial para todos los estados) o Nt=lambda^t*No (notación escalar simple para N total).
#valor reproductivo: sea v*A=lambda*v entonces RV=v1/sum(v1)
M= eigen(t(A)); v= Re(M$vectors[, which.max(Re(M$values))]); RV =v/v[1] #reproductive value v, increases with age.
#sensibilidad y elasticidad: nos dice la importancia relativa de cada transición en determinar lambda. Las elasticidades son la cantidad relativa de los elementos de transición eij=aij*delta_lambda/lambda*delta_aij. con delta_lambda/delta_aij=vij*wij/v*w. Con ambas podemos ver qué controla la tasa de cercimiento de un estado o edad de la población. Si bien no nos dice qué estado y transición está influenciado por el fenómeno, nos permiten predecir los efectos en lambda de un cambio proporcional en la tasa demográdica P o F.
vw.s= v %*% t(w); (S= vw.s/as.numeric(v %*% w)) #calculamos las sensibilidades: la transición más importante es s21.
elas=(A/L1)*S; round(elas,3)  #las elasticidades son sensibilidades sopesadas por las probabilidades de transición, suman cero por lo que nos permiten comparar con otras matrices.
#podríamos discutir más detalles demográficos

#constraste de los modelos con datos reales: derivar intervalos de confianza en los parámetros importantes. Remuestramos los datos de la fila para encontrar los límites de confianza de lambda: bootstrapping (ver Davison & Hinkley). Este método es muy útil en ecología donde los datos no tienen claramente distribuciones paramétricas o modelos demográficos o modelos nulos, donde las aproximaciones analíticas son difíciles o imposibles.
    #bootstrapping: 1) calcular los parámetro(s) observado(s) de interés (e.g. media o lambda) con el modelo y datos, 2) remuestrar nuestras datos con reemplazo para creear un gran número de conjunto de datos y recalcular nuestros parámetros para cada conjunto de datos remuestrado para generar una distribución del parámetro(s) bootstrapped, 3) calcular el intervalo de confianza para el(los) parámetro(s) (intervalo de confianza empírico).
    data(stagedat); data(fruitdat); data(seeddat) #ejemplo: Chamaedorea radicalis.
    str(stagedat); table(stagedat[["Y2004"]]); str(fruitdat); table(fruitdat[["Stage"]], fruitdat[["Y2004"]]); table(seeddat)
    mat1=matrix(0, nrow = 5, ncol = 5); ferts =tapply(fruitdat$Y2004, fruitdat$Stage, mean)/2; mat1[1, 4]=ferts[1]; mat1[1, 5]=ferts[2] #estimando la matriz
    seed.freqs= table(seeddat[, 1]); seedfates= seed.freqs/length(seeddat[, 1]);  mat1[1, 1]=seedfates[2]; mat1[2, 1]= seedfates[3]
    for (i in 2:5) { for (j in 2:5){ mat1[i, j]= {  x= subset(stagedat, stagedat$Y2003 == j); jT= nrow(x); iT=sum(x$Y2004 == i); iT/jT } };  round(mat1, 2)
    ProjMat(stagedat, fruitdat, seeddat) #función para todas las transiciones
    str(DemoInfo(mat1)) # análisis propio
    nL= nrow(stagedat); nF= nrow(fruitdat); nS= nrow(seeddat); n= 5 #boostrapping de la matriz demográfica
    out= lapply(1:n, function(i) { stageR= stagedat[sample(1:nL, nL, replace = TRUE), ]; fruitR= fruitdat[sample(1:nF, nF, replace = TRUE), ]; seedR= as.data.frame(seeddat[sample(1:nS, nS, replace = TRUE), ]); matR= ProjMat(stagedat = stageR, fruitdat = fruitR, seeddat = seedR); DemoInfo(matR) });  sapply(out, function(x) x$lambda) #tenemos 5 estimaciones diferentes de lambda
    estims= DemoInfo(ProjMat(stagedat, fruitdat, seeddat)); estims$lambda; round(estims$Elasticities, 4) #análisis demográfico
    system.time(out.boot= DemoBoot(stagedat, fruitdat, seeddat, n = 1000)); lambdas= sapply(out.boot, function(out) out$lambda); hist(lambdas, prob = T); lines(density(lambdas)) #utilizamos DemoBoot para bootstrap nuestros IC de lambda.
    alpha= 0.05; quantile(lambdas, c(alpha/2, 0.5, 1 - alpha/2))  #para dar un IC real decidimos alfa y calculamos cuantiles. vemos que el IC al 95% (.025 y .975) no incluye el 1. si el ambiente permanece igual la población crecerá lambda>1 con tasa asintótica.
    bias= mean(lambdas) - estims$lambda;  quantile(lambdas - bias, c(alpha/2, 0.5, 1 - alpha/2)) #este boostrap percentil o básico puede dar estimaciones e inferencias inapropiadas. un número de refinamientos pueden hacer para hacer boostrap más preciso y adecuado. Cuando los datos son bastante simétricos, las inferencias pueden ser relativamente válidas, pero cuando existe asimetría puede haber tendencia (la media del boostrap difiere de nuestra estimación observada). Ajustamos el boostrapped para esta tendencia. La tendencia calculada es muy pequeña por lo cual nuestros IC son bastante buenos. Sin embargo, podemos corregir nuestras muestras y restamos la tendencia al lambda boostrapped para obtener el IC. Los cuantiles corregidos indican que la población tiene lambda>1.