A=matrix(c(0,40,70,15,0.3,0,0,0,0,0.4,0,0,0,0,0.55,0), nr=4, byrow = T)
A
## [,1] [,2] [,3] [,4]
## [1,] 0.0 40.0 70.00 15
## [2,] 0.3 0.0 0.00 0
## [3,] 0.0 0.4 0.00 0
## [4,] 0.0 0.0 0.55 0
N0=matrix(c(500,200,350,75), ncol=1)
N0
## [,1]
## [1,] 500
## [2,] 200
## [3,] 350
## [4,] 75
years=15
n.projections=matrix(0,nrow=nrow(A), ncol=years+1)
n.projections
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
n.projections[,1]=N0
n.projections
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 500 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 200 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 350 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 75 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
for (i in 1:years) n.projections[,i+1]=A%*%n.projections[,i]
n.projections
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 500 33625.0 14487.5 408360.00 456795.00 5055303.750 8926106.62
## [2,] 200 150.0 10087.5 4346.25 122508.00 137038.500 1516591.12
## [3,] 350 80.0 60.0 4035.00 1738.50 49003.200 54815.40
## [4,] 75 192.5 44.0 33.00 2219.25 956.175 26951.76
## [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 64904999.40 150030058 858844039 2354399537 11630636907 35615614132
## [2,] 2677831.99 19471500 45009017 257653212 706319861 3489191072
## [3,] 606636.45 1071133 7788600 18003607 103061285 282527944
## [4,] 30148.47 333650 589123 4283730 9901984 56683707
## [,14] [,15] [,16]
## [1,] 160194854594 527415575143 2.233024e+12
## [2,] 10684684240 48058456378 1.582247e+11
## [3,] 1395676429 4273873696 1.922338e+10
## [4,] 155390369 767622036 2.350631e+09
matplot(0:years, t(n.projections), type = "l", lty = 1:4,col = 1, ylab = "Abundancia", xlab = "Año")
legend("topleft", legend = c("Semilla", "Adulto pequeño", "Adulto mediano", "Adulto grande"),lty = 1:4, col = 1, bty = "n")
# Tasa de crecimiento anual (El lambda general de la población es el valor que converge en el tiempo), lambda=Nt+1/Nt para cada año.
# Primero necesitamos la suma de todos los estados aplicando la función sum a cada columna.
n.totals=(apply(n.projections,2,sum))
n.totals
## [1] 1.125000e+03 3.404750e+04 2.467900e+04 4.167742e+05 5.832608e+05
## [6] 5.242302e+06 1.052446e+07 6.821962e+07 1.709063e+08 9.122308e+08
## [11] 2.634340e+09 1.244992e+10 3.944402e+10 1.724306e+11 5.805155e+11
## [16] 2.412822e+12
# Ahora obtenemos cada Lambda dividiendo todos los Nt+1 por cada Nt.
lambdas=n.totals[-1]/n.totals[-(years+1)]
lambdas
## [1] 30.2644444 0.7248403 16.8878095 1.3994645 8.9879211 2.0076038
## [7] 6.4820033 2.5052375 5.3376064 2.8878000 4.7260109 3.1682145
## [13] 4.3715275 3.3666618 4.1563443
#Tenemos un lamba menos que años, ploteamos cada lambda en cada t en vez de en cada t+1
plot(0:(years-1),lambdas,type="b",xlab = "Year",ylab = "L")
#Los lambda que obtuvimos no convergen en el tiempo entonces necesitamos más unidades de tiempo
t=100
Nt=N0/sum(N0)
L.t=numeric(t)
#Creamos una iteración (for-loop) que reusa Nt para cada paso del tiempo, asegurandonos que tengamos un vector numerico vacio para sostener las salidas
for (i in 1:t) L.t[i] <- {
Nt1 <- A %*% Nt
L <- sum(Nt1)/sum(Nt)
Nt <- Nt1/sum(Nt1)
L
}
#Comparamos el resultado directamente al punto estimado de lambda
par(mar = c(5, 4, 3, 2))
plot(1:t, L.t, type = "b", main = quote("Convergencia hacia " * lambda))
L.t
## [1] 30.2644444 0.7248403 16.8878095 1.3994645 8.9879211 2.0076038
## [7] 6.4820033 2.5052375 5.3376064 2.8878000 4.7260109 3.1682145
## [13] 4.3715275 3.3666618 4.1563443 3.5036366 4.0220371 3.5965587
## [19] 3.9367487 3.6588581 3.8819931 3.7002974 3.8465928 3.7277166
## [25] 3.8236021 3.7457958 3.8086270 3.7576892 3.7988542 3.7655014
## [31] 3.7924685 3.7706278 3.7882926 3.7739896 3.7855603 3.7761931
## [37] 3.7837719 3.7776372 3.7826011 3.7785833 3.7818345 3.7792031
## [43] 3.7813325 3.7796091 3.7810038 3.7798750 3.7807885 3.7800492
## [49] 3.7806475 3.7801633 3.7805552 3.7802380 3.7804947 3.7802870
## [55] 3.7804551 3.7803190 3.7804291 3.7803400 3.7804121 3.7803538
## [61] 3.7804010 3.7803628 3.7803937 3.7803687 3.7803889 3.7803725
## [67] 3.7803858 3.7803751 3.7803838 3.7803767 3.7803824 3.7803778
## [73] 3.7803815 3.7803785 3.7803810 3.7803790 3.7803806 3.7803793
## [79] 3.7803803 3.7803795 3.7803802 3.7803796 3.7803801 3.7803797
## [85] 3.7803800 3.7803798 3.7803800 3.7803798 3.7803799 3.7803798
## [91] 3.7803799 3.7803798 3.7803799 3.7803799 3.7803799 3.7803799
## [97] 3.7803799 3.7803799 3.7803799 3.7803799
#Usando ánalysis de valores y vectores "propios" (Eigenanalysis) de la matriz llegamos directamente a la convergencia de lambda
eigs.A=eigen(A)
eigs.A
## $values
## [1] 3.7803799 -3.0594388 -0.5710462 -0.1498949
##
## $vectors
## [,1] [,2] [,3] [,4]
## [1,] 0.996830360 -0.995143171 -0.8066003 -0.04893724
## [2,] 0.079105571 0.097580952 0.4237487 0.09794314
## [3,] 0.008370119 -0.012758020 -0.2968227 -0.26136491
## [4,] 0.001217752 0.002293529 0.2858832 0.95901026
dom.pos=which.max(eigs.A[["values"]])
L1 <- Re(eigs.A[["values"]][dom.pos])
L1
## [1] 3.78038
points(t, L1, pch = 19, cex = 1.5)
#La distribución estable de estados se alcanza cuando la abundancia relativa de los estados es la misma
w <- Re(eigs.A[["vectors"]][, dom.pos])
ssd <- w/sum(w)
round(ssd, 3)
## [1] 0.918 0.073 0.008 0.001
# Lo anterior muestra que si la matriz de projección no cambia en el tiempo, estara compuesta de los porcentajes calculados para cada estado
#El valor reproductivo es el número de individuos a nacer de una determinada edad, lo calculamos con la matriz transpuesta de la matriz de proyecciones (A)
M <- eigen(t(A))
v <- Re(M$vectors[, which.max(Re(M$values))])
v
## [1] 0.0430286 0.5422149 0.8215853 0.1707313
RV <- v/v[1]
round(RV,3)
## [1] 1.000 12.601 19.094 3.968
#La importancia relativa de cada transición del ciclo de vida se hace calculando la "sensibilidad" y "elasticidad" de la matriz
#Las sensibilidades son la contribución directa de cada transición a lambda
vw.s <- v %*% t(w)
(S <- vw.s/as.numeric(v %*% w))
## [,1] [,2] [,3] [,4]
## [1,] 0.4618567 0.03665161 0.003878087 0.0005642153
## [2,] 5.8199788 0.46185667 0.048868810 0.0071098267
## [3,] 8.8186603 0.69982335 0.074047939 0.0107730884
## [4,] 1.8325804 0.14542827 0.015387689 0.0022387245
#Las elasticidades son sensibilidades ponderadas por las probabilidades de transición
elas <- (A/L1) * S
round(elas, 3)
## [,1] [,2] [,3] [,4]
## [1,] 0.000 0.388 0.072 0.002
## [2,] 0.462 0.000 0.000 0.000
## [3,] 0.000 0.074 0.000 0.000
## [4,] 0.000 0.000 0.002 0.000
#Dos propiedades importantes tienen estas elasticidades, transiciones imposibles tienen elasticidades iguales a cero y todas las elasticidades de la matriz suman uno
#Así es facil comparar elasticidades para diferentes matrices y diferentes organismos
#Lo anterior es importante en el manejo de especies invasoras o en peligro donde buscamos tener el maximo impacto con el mínimo esfuerzo