Project
Colin
April 4, 2019
Installing and loading the packages
library(popbio)Building the initial Matrix
s0<-0.0384
s1<-0.83
s2<-0.83
s3<-0.83
s4<-0.83
s5<-0.83
s6<-0.83
s7<-0.83
s8<-0.83
s9<-0.83
s10<-0.83
s11<-0.83
s12<-0.83
s13<-0.83
s14<-0.83
s15<-0.7425
s16<-0.7425
s17<-0.7425
s18<-0.7425
s19<-0.7425
s20<-0.86
f0<-0
f1<-0
f2<-0
f3<-41.61356
projection_matrix<-matrix(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,41.61356,
s0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,s1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,s2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,s3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,s4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,s5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,s6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,s7,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,s8,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,s9,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,s10,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,s11,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,s12,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,s13,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,s14,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,s15,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,s16,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,s17,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,s18,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,s19,s20),
nrow=21, byrow=T)
lambda(projection_matrix)## [1] 0.9440294Building the Stage classified model
Now we need to convert all this to the stage-classified model (see word document for more detailed explanation)
Juvenile to sub-adult transition
y here represents the annual probability of an individual growing to the next STAGE (i.e. from juvenile to sub-adult)
surv_juv<-0.83
tort_lam<-lambda(projection_matrix)
y_juv<- (((surv_juv/tort_lam)^14) - ((surv_juv/tort_lam)^13)) / (((surv_juv/tort_lam)^14) - 1)
y_juv ## [1] 0.02713383So, now we have to calculate the probability of of surviving and remaining in each stage
P_juv<- surv_juv * (1 - y_juv)
P_juv## [1] 0.8074789We also need to calculate the probability of surviving and growing to stage i+1
G_juv<- surv_juv * y_juv
G_juv## [1] 0.02252108Sub-adult to adult transition
Now we can do the same thing, except with the transition from sub-adult to adult status
surv_subadult<-0.7425
tort_lam<-lambda(projection_matrix) #lambda is obviously still the same here
y_subadult<-(((surv_subadult/tort_lam)^5) - ((surv_subadult/tort_lam)^4)) / (((surv_subadult/tort_lam)^5) - 1)
y_subadult ## [1] 0.1168731P_subadult<- surv_subadult * (1 - y_subadult)
P_subadult## [1] 0.6557218G_subadult<- surv_subadult * y_subadult
G_subadult## [1] 0.08677824Building more intuitive matrix
Now, we can set up the more intuitive matrix
f2<- 0 #reproductive output at stage 2
f3<- 0 #reproductive output at stage 3
f4<- 41.61356 #reproductive output at stage 4
g1<- 0.0384 # probability of surviving and transitioning to the next stage (here from hatchlin to juvenile)
g2<- 0.02252108
g3<- 0.08677824
p1<- 0 #probability of surviving and staying at age 1 (this is 0 because all individuals that survive transition to juv status)
p2<- 0.8074789
p3<- 0.6557218
p4<- 0.86
turt_matrix<- matrix(c(p1,f2,f3,f4,
g1,p2,0,0,
0,g2,p3,0,
0,0,g3,p4),
nrow = 4, byrow = T)
lambda(turt_matrix) #As you can see, the lambda remains unchanged from what it was## [1] 0.9440294Sensitivity / elasticity analysis
First I am just setting up the matrix of vital rates with the names of the stage states, and then checking out the elasticity/ sensitivity of each of them
rownames(turt_matrix)<-c('egg/hatchling', 'juvenile', 'sub-adult', 'adult')
colnames(turt_matrix)<-c('egg/hatchling', 'juvenile', 'sub-adult', 'adult')
elasticity(turt_matrix)## egg/hatchling juvenile sub-adult adult
## egg/hatchling 0.00000000 0.00000000 0.00000000 0.04459845
## juvenile 0.04459845 0.26372885 0.00000000 0.00000000
## sub-adult 0.00000000 0.04459845 0.10143394 0.00000000
## adult 0.00000000 0.00000000 0.04459845 0.45644342sensitivity(turt_matrix)## egg/hatchling juvenile sub-adult adult
## egg/hatchling 0.04459845 0.01254174 0.0009796948 0.001011743
## juvenile 1.09641266 0.30832730 0.0240849143 0.024872801
## sub-adult 6.64780265 1.86945947 0.1460323857 0.150809523
## adult 22.08632041 6.21099679 0.4851705488 0.501041866Barplot sensitivity and elasticity
# These commands just get the vital rates ready to be stored in a way that works with the 'vitalsense' command
tort.el<-expression(0,0,0,f4,
g1,p2,0,0,
0,g2,p3,0,
0,0,g3,p4)
tort.vr<-list(f4=41.61356, g1=0.0384, p2=0.8074789,
g2=0.02252108, p3=0.6557218, g3=0.08677824, p4=0.86)
# This creates a matrix with the vital rates and their corresponding elasticity and sensitivity
x<-vitalsens(tort.el, tort.vr)
# This is just the actual barplot
barplot(t(x[,2:3]), beside=TRUE, legend=TRUE, las=1, xlab="Vital rate",
main="Caretta caretta vital rate sensitivity and elasticity")
abline(h=0)
Plotting the initial Decline
Here, I am plotting how long it will take each stage class to go down to 0, based upon the current value of lambda (I copied and pasted all this code from a website; I just changed the data) Remember that you can change the matrix values to simulate what might happen if you reduce hatchling mortality, or increase female fecundity, etc. etc.
stages<-c("egg/hatchling", "juvenile", "sub-adult","adult")
projection_matrix<-matrix(c(0,0,0,41.613,
0.0384,0.8074789,0,0,
0,0.02252108,0.6557218,0,
0,0,0.08677824,0.86),
nrow=4, byrow=T,dimnames=list(stages,stages))
Abundance_year0<-c(576.5,2813.5,11254,2813.5) # I used these values just to make it intuitive; you can use whatever values you want
nYears <- 50 # set the number of years to project
TMat<-projection_matrix # define the projection matrix
InitAbund <- Abundance_year0 # define the initial abundance
Year1 <- projection_matrix %*% Abundance_year0
## NOTE: the code below can be re-used without modification:
allYears <- matrix(0,nrow=nrow(TMat),ncol=nYears+1) # build a storage array for all abundances!
allYears[,1] <- InitAbund # set the year 0 abundance
for(t in 2:(nYears+1)){ # loop through all years
allYears[,t] <- TMat %*% allYears[,t-1]
}
plot(1,1,pch="",ylim=c(0,max(allYears)),xlim=c(0,nYears+1),xlab="Years",ylab="Abundance",xaxt="n") # set up blank plot
cols <- rainbow(ncol(TMat)) # set up colors to use
for(s in 1:ncol(TMat)){
points(allYears[s,],col=cols[s],type="l",lwd=2) # plot out each life stage abundance, one at a time
}
axis(1,at=seq(1,nYears+1),labels = seq(0,nYears)) # label the axis
legend("topleft",col=cols,lwd=rep(2,ncol(TMat)),legend=paste("Stage ",seq(1:ncol(TMat)))) # put a legend on the plot
Stable Age structure
Now to graph how long it takes to reach stable age structure
stages<-c("egg/hatchling", "juvenile", "sub-adult","adult")
projection_matrix<-matrix(c(0,0,0,41.613,
0.0384,0.8074789,0,0,
0,0.02252108,0.6557218,0,
0,0,0.08677824,0.86),
nrow=4, byrow=T,dimnames=list(stages,stages))
dimnames=list(stages,stages)
n<-c(576.5,0,0,0) # These are the initial abundance values, used just as a better illustration
p<-pop.projection(projection_matrix,n, 30) ‘pop.projection’ asks for the projection matrix, the vector of initial stage abundances, and the number of years to simulate
stage.vector.plot(p$stage.vectors, col=2:4) #plots the simulation
Here, remember in class how ALL populations eventually hit a stable age structure (as shown above), but we can mess with the initial abundance values to change how long it takes to reach stable age structure.