setwd("~/Dropbox/Clases_UPJ/Biología de poblaciones/Alumnos_BioP/2018_10/Clases/Prácticas/Práctica 2")
# Sección I
# Punto 1a, cálculo de r
# Cargamos la base de datos
datos1<-read.table("invasive.csv", sep = ",", header=T)
attach(datos1)
# Transformamos los datos del numero de individuos a logaritmo natural y el tiempo en fechas lo convertimos a número de años transcurridos
datos1.tra<-transform(datos1,n=log(n), tiempo=tiempo-1962)
detach(datos1)
attach(datos1.tra)
# Ajustamos un modelo lineal al número de individuos en función del tiempo
fit1<-lm(n~tiempo)
summary(fit1)
##
## Call:
## lm(formula = n ~ tiempo)
##
## Residuals:
## 1 2 3 4 5 6
## -0.69394 1.46597 -0.32546 -0.67268 0.29003 -0.06393
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.38709 0.75809 1.830 0.14127
## tiempo 0.18810 0.02986 6.299 0.00325 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9052 on 4 degrees of freedom
## Multiple R-squared: 0.9084, Adjusted R-squared: 0.8855
## F-statistic: 39.68 on 1 and 4 DF, p-value: 0.003246
# Punto 1b, cálculo para el tiempo doble usando el r estimado
log(2)/0.18810
## [1] 3.684993
# Punto 1c,despejamos el tiempo en la ecuación lineal y creamos una función para predecir cuanto demorará la población en alcanzar 10000 individuos
m.time1=function(r,N=10000){(log(N)-1.38709)/r}
#Creamos un vector con diferentes valores para r que nos permita predecir en cuanto tiempo la población alcanza 10000 individuos si crece 2, 3 o 4 veces mas rápido
rs=c(0.18810,0.18810*2,0.18810*3,0.18810*4)
rs
## [1] 0.1881 0.3762 0.5643 0.7524
m.time1(rs)
## [1] 41.59091 20.79546 13.86364 10.39773
detach(datos1.tra)
#Sección II
#Punto 5a, Matriz de Leslie
A=matrix(c(0,40,50,20,0.25,0,0,0,0,0.4,0,0,0,0,0.25,0), nr=4, byrow = T)
A
## [,1] [,2] [,3] [,4]
## [1,] 0.00 40.0 50.00 20
## [2,] 0.25 0.0 0.00 0
## [3,] 0.00 0.4 0.00 0
## [4,] 0.00 0.0 0.25 0
N0=matrix(c(400,100,250,75), ncol=1)
N0
## [,1]
## [1,] 400
## [2,] 100
## [3,] 250
## [4,] 75
# Punto 5b, proyección a 10 años
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
n.projections[,1]=N0
n.projections
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 400 0 0 0 0 0 0 0 0 0 0
## [2,] 100 0 0 0 0 0 0 0 0 0 0
## [3,] 250 0 0 0 0 0 0 0 0 0 0
## [4,] 75 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] [,8]
## [1,] 400 18000.0 7250 182200.0 162700 1867250.00 2541625.0 19577100.0
## [2,] 100 100.0 4500 1812.5 45550 40675.00 466812.5 635406.2
## [3,] 250 40.0 40 1800.0 725 18220.00 16270.0 186725.0
## [4,] 75 62.5 10 10.0 450 181.25 4555.0 4067.5
## [,9] [,10] [,11]
## [1,] 34833850.00 209412750.00 447494812.5
## [2,] 4894275.00 8708462.50 52353187.5
## [3,] 254162.50 1957710.00 3483385.0
## [4,] 46681.25 63540.62 489427.5
#Punto 5c, gráfica por estados
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")
#Punto 6a, cálculo de elasticidades
#Usando ánalisis de valores y vectores "propios" (Eigenanalysis) de la matriz llegamos directamente a la convergencia de lambda sin tenerlo que hacer paso a paso.
eigs.A=eigen(A)
eigs.A
## $values
## [1] 3.3936354 -2.8857187 -0.3698824 -0.1380344
##
## $vectors
## [,1] [,2] [,3] [,4]
## [1,] 0.9972599647 -0.99619648 -0.6688791 -0.0904684
## [2,] 0.0734654616 0.08630402 0.4520891 0.1638512
## [3,] 0.0086592050 -0.01196292 -0.4889004 -0.4748129
## [4,] 0.0006379004 0.00103639 0.3304432 0.8599541
dom.pos=which.max(eigs.A[["values"]])
L1 <- Re(eigs.A[["values"]][dom.pos])
L1
## [1] 3.393635
#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.923 0.068 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.04713749 0.63986977 0.71496312 0.27779935
RV <- v/v[1]
round(RV,3)
## [1] 1.000 13.575 15.168 5.893
#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.468281 0.0344970 0.004066082 0.0002995373
## [2,] 6.356700 0.4682810 0.055195202 0.0040660822
## [3,] 7.102704 0.5232371 0.061672758 0.0045432664
## [4,] 2.759760 0.2033041 0.023962988 0.0017652889
#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.407 0.060 0.002
## [2,] 0.468 0.000 0.000 0.000
## [3,] 0.000 0.062 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