A=matrix(c(0,60,20,5,0.3,0,0,0,0,0.4,0,0,0,0,0.25,0), nr=4, byrow = T)
A
## [,1] [,2] [,3] [,4]
## [1,] 0.0 60.0 20.00 5
## [2,] 0.3 0.0 0.00 0
## [3,] 0.0 0.4 0.00 0
## [4,] 0.0 0.0 0.25 0
N0=matrix(c(400,100,200,50), ncol=1)
N0
## [,1]
## [1,] 400
## [2,] 100
## [3,] 200
## [4,] 50
years=20
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] [,17] [,18] [,19] [,20] [,21]
## [1,] 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0
## [4,] 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] [,12] [,13]
## [1,] 400 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 100 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 200 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 50 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0
## [4,] 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] [,8]
## [1,] 400 10250 8250 185510 173160.0 3360517.5 3563341.5 60932725.5
## [2,] 100 120 3075 2475 55653.0 51948.0 1008155.2 1069002.4
## [3,] 200 40 48 1230 990.0 22261.2 20779.2 403262.1
## [4,] 50 50 10 12 307.5 247.5 5565.3 5194.8
## [,9] [,10] [,11] [,12] [,13]
## [1,] 72231363.0 1105845156.2 1446937576 20087707992 28709739455
## [2,] 18279817.6 21669408.9 331753547 434081273 6026312398
## [3,] 427601.0 7311927.1 8667764 132701419 173632509
## [4,] 100815.5 106900.2 1827982 2166941 33175355
## [,14] [,15] [,16] [,17] [,18]
## [1,] 365217270814 565202850008 6.645827e+12 1.105448e+13 1.210362e+14
## [2,] 8612921837 109565181244 1.695609e+11 1.993748e+12 3.316344e+12
## [3,] 2410524959 3445168735 4.382607e+10 6.782434e+10 7.974993e+11
## [4,] 43408127 602631240 8.612922e+08 1.095652e+10 1.695609e+10
## [,19] [,20] [,21]
## [1,] 2.150154e+14 2.206179e+15 4.162422e+15
## [2,] 3.631085e+13 6.450462e+13 6.618536e+14
## [3,] 1.326538e+12 1.452434e+13 2.580185e+13
## [4,] 1.993748e+11 3.316344e+11 3.631085e+12
# Gráfica por estados
matplot(0:years, t(n.projections), type = "l", lty = 1:3,col = 1, ylab = "Abundancia", xlab = "Año")
legend("topleft", legend = c("Semillas", "Adultos pequeños", "Adultos medianos", "Adultos grandes"),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] 7.500000e+02 1.046000e+04 1.138300e+04 1.892270e+05 2.301105e+05
## [6] 3.434974e+06 4.597841e+06 6.241018e+07 9.103960e+07 1.134933e+09
## [11] 1.789187e+09 2.065666e+10 3.494286e+10 3.762841e+11 6.788158e+11
## [16] 6.860076e+12 1.312701e+13 1.251670e+14 2.528522e+14 2.285539e+15
## [21] 4.853709e+15
# Ahora obtenemos cada Lambda dividiendo todos los Nt+1 por cada Nt.
lambdas=n.totals[-1]/n.totals[-(years+1)]
lambdas
## [1] 13.946667 1.088241 16.623649 1.216055 14.927499 1.338537 13.573802
## [8] 1.458730 12.466371 1.576469 11.545277 1.691603 10.768556 1.803998
## [15] 10.105945 1.913537 9.535071 2.020119 9.039034 2.123660
#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, la covergencia se da en el año 606
t=700
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] 13.946667 1.088241 16.623649 1.216055 14.927499 1.338537 13.573802
## [8] 1.458730 12.466371 1.576469 11.545277 1.691603 10.768556 1.803998
## [15] 10.105945 1.913537 9.535071 2.020119 9.039034 2.123660 8.604828
## [22] 2.224093 8.222282 2.321366 7.883328 2.415443 7.581480 2.506303
## [29] 7.311468 2.593936 7.068963 2.678348 6.850376 2.759556 6.652708
## [36] 2.837587 6.473431 2.912479 6.310401 2.984279 6.161785 3.053040
## [43] 6.026011 3.118825 5.901718 3.181701 5.787725 3.241740 5.683002
## [50] 3.299020 5.586645 3.353621 5.497859 3.405626 5.415941 3.455120
## [57] 5.340268 3.502191 5.270286 3.546925 5.205498 3.589410 5.145463
## [64] 3.629734 5.089781 3.667984 5.038095 3.704246 4.990080 3.738606
## [71] 4.945445 3.771145 4.903924 3.801946 4.865275 3.831089 4.829280
## [78] 3.858650 4.795739 3.884705 4.764468 3.909326 4.735300 3.932585
## [85] 4.708083 3.954548 4.682675 3.975282 4.658947 3.994849 4.636780
## [92] 4.013309 4.616065 4.030721 4.596702 4.047139 4.578595 4.062616
## [99] 4.561660 4.077204 4.545816 4.090949 4.530990 4.103899 4.517113
## [106] 4.116097 4.504122 4.127584 4.491958 4.138401 4.480566 4.148583
## [113] 4.469896 4.158168 4.459900 4.167190 4.450534 4.175679 4.441757
## [120] 4.183666 4.433531 4.191181 4.425821 4.198251 4.418593 4.204900
## [127] 4.411816 4.211154 4.405462 4.217035 4.399504 4.222565 4.393917
## [134] 4.227765 4.388677 4.232654 4.383761 4.237250 4.379151 4.241571
## [141] 4.374826 4.245633 4.370768 4.249450 4.366962 4.253038 4.363390
## [148] 4.256411 4.360038 4.259580 4.356894 4.262558 4.353943 4.265357
## [155] 4.351173 4.267987 4.348574 4.270458 4.346135 4.272780 4.343846
## [162] 4.274961 4.341697 4.277011 4.339681 4.278936 4.337788 4.280745
## [169] 4.336011 4.282445 4.334343 4.284041 4.332777 4.285541 4.331308
## [176] 4.286950 4.329928 4.288274 4.328633 4.289517 4.327417 4.290685
## [183] 4.326275 4.291782 4.325204 4.292812 4.324198 4.293780 4.323253
## [190] 4.294689 4.322366 4.295543 4.321533 4.296345 4.320752 4.297098
## [197] 4.320018 4.297805 4.319329 4.298470 4.318682 4.299094 4.318074
## [204] 4.299680 4.317504 4.300230 4.316968 4.300747 4.316465 4.301233
## [211] 4.315993 4.301689 4.315549 4.302117 4.315133 4.302519 4.314742
## [218] 4.302897 4.314375 4.303252 4.314030 4.303585 4.313707 4.303898
## [225] 4.313403 4.304192 4.313118 4.304468 4.312850 4.304727 4.312598
## [232] 4.304970 4.312362 4.305199 4.312140 4.305413 4.311932 4.305615
## [239] 4.311736 4.305804 4.311552 4.305982 4.311380 4.306149 4.311218
## [246] 4.306306 4.311066 4.306453 4.310923 4.306591 4.310789 4.306721
## [253] 4.310663 4.306843 4.310545 4.306958 4.310434 4.307065 4.310329
## [260] 4.307166 4.310231 4.307261 4.310139 4.307350 4.310053 4.307434
## [267] 4.309972 4.307512 4.309896 4.307586 4.309824 4.307655 4.309757
## [274] 4.307720 4.309694 4.307781 4.309635 4.307839 4.309579 4.307892
## [281] 4.309527 4.307943 4.309478 4.307990 4.309432 4.308035 4.309389
## [288] 4.308077 4.309348 4.308116 4.309310 4.308153 4.309274 4.308188
## [295] 4.309241 4.308220 4.309209 4.308251 4.309179 4.308280 4.309152
## [302] 4.308307 4.309125 4.308332 4.309101 4.308356 4.309078 4.308378
## [309] 4.309056 4.308399 4.309036 4.308419 4.309017 4.308437 4.308999
## [316] 4.308455 4.308982 4.308471 4.308966 4.308486 4.308951 4.308501
## [323] 4.308937 4.308514 4.308924 4.308527 4.308912 4.308539 4.308900
## [330] 4.308550 4.308890 4.308561 4.308879 4.308570 4.308870 4.308580
## [337] 4.308861 4.308588 4.308852 4.308597 4.308844 4.308604 4.308837
## [344] 4.308611 4.308830 4.308618 4.308824 4.308625 4.308817 4.308631
## [351] 4.308812 4.308636 4.308806 4.308641 4.308801 4.308646 4.308796
## [358] 4.308651 4.308792 4.308655 4.308788 4.308659 4.308784 4.308663
## [365] 4.308780 4.308667 4.308776 4.308670 4.308773 4.308673 4.308770
## [372] 4.308676 4.308767 4.308679 4.308764 4.308682 4.308762 4.308684
## [379] 4.308759 4.308687 4.308757 4.308689 4.308755 4.308691 4.308753
## [386] 4.308693 4.308751 4.308695 4.308749 4.308696 4.308748 4.308698
## [393] 4.308746 4.308699 4.308745 4.308701 4.308743 4.308702 4.308742
## [400] 4.308703 4.308741 4.308705 4.308740 4.308706 4.308739 4.308707
## [407] 4.308738 4.308708 4.308737 4.308709 4.308736 4.308709 4.308735
## [414] 4.308710 4.308734 4.308711 4.308734 4.308712 4.308733 4.308712
## [421] 4.308732 4.308713 4.308732 4.308713 4.308731 4.308714 4.308731
## [428] 4.308715 4.308730 4.308715 4.308730 4.308715 4.308729 4.308716
## [435] 4.308729 4.308716 4.308728 4.308717 4.308728 4.308717 4.308728
## [442] 4.308717 4.308727 4.308718 4.308727 4.308718 4.308727 4.308718
## [449] 4.308727 4.308719 4.308726 4.308719 4.308726 4.308719 4.308726
## [456] 4.308719 4.308726 4.308719 4.308725 4.308720 4.308725 4.308720
## [463] 4.308725 4.308720 4.308725 4.308720 4.308725 4.308720 4.308725
## [470] 4.308720 4.308725 4.308720 4.308724 4.308721 4.308724 4.308721
## [477] 4.308724 4.308721 4.308724 4.308721 4.308724 4.308721 4.308724
## [484] 4.308721 4.308724 4.308721 4.308724 4.308721 4.308724 4.308721
## [491] 4.308724 4.308721 4.308723 4.308721 4.308723 4.308722 4.308723
## [498] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [505] 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723
## [512] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [519] 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723
## [526] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [533] 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723
## [540] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [547] 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723
## [554] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [561] 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723
## [568] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [575] 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723
## [582] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [589] 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723
## [596] 4.308722 4.308723 4.308722 4.308723 4.308722 4.308723 4.308722
## [603] 4.308723 4.308722 4.308723 4.308722 4.308722 4.308722 4.308722
## [610] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [617] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [624] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [631] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [638] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [645] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [652] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [659] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [666] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [673] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [680] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [687] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
## [694] 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722 4.308722
#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
## eigen() decomposition
## $values
## [1] 4.30872247+0.00000000i -4.17538096+0.00000000i -0.06667076+0.06239164i
## [4] -0.06667076-0.06239164i
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.9975640650+0i 0.9974051688+0i 0.01504872+0.01838392i
## [2,] -0.0694566015+0i -0.0716632934+0i 0.00517024-0.07788416i
## [3,] -0.0064479996+0i 0.0068653178+0i -0.24966217+0.23363815i
## [4,] -0.0003741248+0i -0.0004110594+0i 0.93617568+0.00000000i
## [,4]
## [1,] 0.01504872-0.01838392i
## [2,] 0.00517024+0.07788416i
## [3,] -0.24966217-0.23363815i
## [4,] 0.93617568+0.00000000i
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] 4.308722
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.929 0.065 0.006 0.000
# 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.06582356 0.94538480 0.30996825 0.07638408
RV <- v/v[1]
round(RV,3)
## [1] 1.000 14.362 4.709 1.160
#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.4923990 0.03428387 0.003182741 0.0001846685
## [2,] 7.0720349 0.49239896 0.045711829 0.0026522844
## [3,] 2.3187450 0.16144542 0.014987776 0.0008696183
## [4,] 0.5713979 0.03978427 0.003693370 0.0002142961
#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.477 0.015 0
## [2,] 0.492 0.000 0.000 0
## [3,] 0.000 0.015 0.000 0
## [4,] 0.000 0.000 0.000 0
#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