A=matrix(c(0,0,0,0,0,322.380,0.966,0,0,0,0,0,0.013,0.010,0.125,0,0,3.448,0.007,0,0.125,0.238,0,30.170, 0.008,0,0,0.245,0.167,0.862,0,0,0,0.023,0.750,0), nr=6, byrow = T)
A
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.000 0.00 0.000 0.000 0.000 322.380
## [2,] 0.966 0.00 0.000 0.000 0.000 0.000
## [3,] 0.013 0.01 0.125 0.000 0.000 3.448
## [4,] 0.007 0.00 0.125 0.238 0.000 30.170
## [5,] 0.008 0.00 0.000 0.245 0.167 0.862
## [6,] 0.000 0.00 0.000 0.023 0.750 0.000
N0=matrix(c(200,150,140,130,110,100), ncol=1)
N0
## [,1]
## [1,] 200
## [2,] 150
## [3,] 140
## [4,] 130
## [5,] 110
## [6,] 100
years=10
n.projections=matrix(0,nrow=nrow(A), ncol=years+1)
n.projections
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0
n.projections[,1]=N0
n.projections
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 200 0 0 0 0 0 0 0 0 0 0
## [2,] 150 0 0 0 0 0 0 0 0 0 0
## [3,] 140 0 0 0 0 0 0 0 0 0 0
## [4,] 130 0 0 0 0 0 0 0 0 0 0
## [5,] 110 0 0 0 0 0 0 0 0 0 0
## [6,] 100 0 0 0 0 0 0 0 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,] 200 32238.00 27560.2662 56110.9869 293968.685 393740.706 954483.146
## [2,] 150 193.20 31141.9080 26623.2171 54203.213 283973.750 380353.522
## [3,] 140 366.40 761.5955 1365.0344 4310.432 9113.664 19306.204
## [4,] 130 3066.84 3580.6072 6391.4643 29595.696 46488.667 104285.201
## [5,] 110 138.02 1106.0215 1432.4696 3040.051 11163.192 18956.060
## [6,] 100 85.49 174.0523 911.8701 1221.356 2960.739 9441.634
## [,8] [,9] [,10] [,11]
## [1,] 3043793.82 5356538.50 13120622.7 34514884.4
## [2,] 922030.72 2940304.83 5174416.2 12674521.6
## [3,] 51179.84 112477.71 253428.8 623143.2
## [4,] 318768.62 604863.75 1423409.3 3692377.1
## [5,] 44490.09 124201.16 246868.3 587215.3
## [6,] 16615.60 40699.25 107062.7 217889.6
# 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] 830.00 36087.95 64324.45 92835.04 386339.43
## [6] 747440.72 1486825.77 4396878.69 9179085.20 20325807.99
## [11] 52310031.19
# Ahora obtenemos cada Lambda dividiendo todos los Nt+1 por cada Nt.
lambdas=n.totals[-1]/n.totals[-(years+1)]
lambdas
## [1] 43.479458 1.782436 1.443231 4.161569 1.934674 1.989222 2.957225
## [8] 2.087637 2.214361 2.573577
#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=60
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] 43.479458 1.782436 1.443231 4.161569 1.934674 1.989222 2.957225
## [8] 2.087637 2.214361 2.573577 2.193714 2.294100 2.424841 2.256220
## [15] 2.318653 2.363861 2.289720 2.324479 2.338667 2.306639 2.324785
## [22] 2.328391 2.314845 2.323915 2.324300 2.318707 2.323095 2.322725
## [29] 2.320479 2.322546 2.322146 2.321275 2.322225 2.321948 2.321624
## [36] 2.322052 2.321887 2.321774 2.321963 2.321873 2.321838 2.321919
## [43] 2.321873 2.321864 2.321898 2.321875 2.321874 2.321888 2.321877
## [50] 2.321878 2.321884 2.321878 2.321879 2.321882 2.321879 2.321880
## [57] 2.321881 2.321880 2.321880 2.321880
#Usando ánalysis de valores y vectores "propios" (Eigenanalysis) de la matriz tambien podemos llegar a la convergencia de lambda
eigs.A=eigen(A)
eigs.A
## $values
## [1] 2.3218801+0.000000i -0.9573856+1.488566i -0.9573856-1.488566i
## [4] 0.1423853+0.198015i 0.1423853-0.198015i -0.1618797+0.000000i
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.918358907+0i 0.875276087+0.000000000i 0.875276087+0.000000000i
## [2,] -0.382076011+0i -0.258422163-0.401800898i -0.258422163+0.401800898i
## [3,] -0.017554688+0i 0.004411771-0.003094969i 0.004411771+0.003094969i
## [4,] -0.099898451+0i 0.073170651-0.010562731i 0.073170651+0.010562731i
## [5,] -0.017413276+0i -0.006947229-0.009994197i -0.006947229+0.009994197i
## [6,] -0.006614304+0i -0.002599345+0.004041522i -0.002599345-0.004041522i
## [,4] [,5]
## [1,] 1.426564e-01+1.983921e-01i 1.426564e-01-1.983921e-01i
## [2,] 9.678394e-01+0.000000e+00i 9.678394e-01+0.000000e+00i
## [3,] 2.093852e-02-5.537957e-02i 2.093852e-02+5.537957e-02i
## [4,] -4.648716e-03-7.049398e-03i -4.648716e-03+7.049398e-03i
## [5,] 8.511912e-05+2.339139e-04i 8.511912e-05-2.339139e-04i
## [6,] -5.885111e-05+1.752473e-04i -5.885111e-05-1.752473e-04i
## [,6]
## [1,] -1.652062e-01+0i
## [2,] 9.858505e-01+0i
## [3,] -2.787531e-02+0i
## [4,] 5.346755e-03+0i
## [5,] -1.818725e-04+0i
## [6,] 8.295651e-05+0i
dom.pos=which.max(eigs.A[["values"]])
## Warning in which.max(eigs.A[["values"]]): imaginary parts discarded in
## coercion
L1 <- Re(eigs.A[["values"]][dom.pos])
L1
## [1] 2.32188
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.637 0.265 0.012 0.069 0.012 0.005
# 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] 1.299554e-03 1.201025e-05 2.788635e-03 4.901038e-02 3.283103e-01
## [6] 9.432926e-01
RV <- v/v[1]
round(RV,3)
## [1] 1.000 0.009 2.146 37.713 252.633 725.858
#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] [,5]
## [1,] 6.593963e-02 0.027433664 1.260454e-03 7.172867e-03 1.250301e-03
## [2,] 6.094021e-04 0.000253537 1.164889e-05 6.629034e-05 1.155505e-05
## [3,] 1.414959e-01 0.058868245 2.704733e-03 1.539182e-02 2.682945e-03
## [4,] 2.486796e+00 1.034611834 4.753580e-02 2.705119e-01 4.715287e-02
## [5,] 1.665853e+01 6.930649121 3.184324e-01 1.812103e+00 3.158673e-01
## [6,] 4.786284e+01 19.912957587 9.149116e-01 5.206487e+00 9.075415e-01
## [,6]
## [1,] 4.749175e-04
## [2,] 4.389102e-06
## [3,] 1.019097e-03
## [4,] 1.791067e-02
## [5,] 1.199798e-01
## [6,] 3.447229e-01
#Las elasticidades son sensibilidades ponderadas por las probabilidades de transición
elas <- (A/L1) * S
round(elas, 3)
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.000 0 0.000 0.000 0.000 0.066
## [2,] 0.000 0 0.000 0.000 0.000 0.000
## [3,] 0.001 0 0.000 0.000 0.000 0.002
## [4,] 0.007 0 0.003 0.028 0.000 0.233
## [5,] 0.057 0 0.000 0.191 0.023 0.045
## [6,] 0.000 0 0.000 0.052 0.293 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