# sensitivity and elasticity to lower-level vital rates
library(popbio)
?popbio # to see the menu
# Data from:
# Hiraldo, F., J.J. Negro, J.A. Donazar, & P. Gaona. 1996.
# A demographic model for a population of the endangered
# lesser kestrel in southern Spain.
# Journal of Applied Ecology 33(5):1085-1093.
## Based on Table 6 - See Figure 1 for definitions of rates
kest.vr<- list(b = 0.9321, c0 = 0.3847, ca = 0.925, s0 = 0.3409, sa = 0.7101)
kest.el <- expression( c0*b*s0, ca*b*s0, sa, sa)
# store a version of the matrix as symbols so can check it
A.sym<-matrix(as.character(kest.el), nrow=sqrt(length(kest.el)), byrow=TRUE)
A.sym
# use the rates and expressions to create the matrix with numbers
A<-matrix(sapply(kest.el, eval,kest.vr, NULL),
nrow=sqrt(length(kest.el)), byrow=TRUE)
A
# get the s & e values associated with each matrix element
E=eigen.analysis(A)
E
# check that the reproductive values are the same as those
# reported in the paper (the values in the paper have
# obviously been rounded)
.69/.69 # RV for juveniles = 1
.82/.69 # RV for adults relative to that of juveniles
# get the s & e values associated with vital rates
x<-vitalsens(kest.el, kest.vr)
x
# compare these to what's reported in Table 7
# Which is correct: R or the paper?
# Let's change 1 rate 'by hand' and see what we get
kest.vr2<- list(b = 0.9321+.001, c0 = 0.3847, ca = 0.925, s0 = 0.3409, sa = 0.7101)
A.mod.b<-matrix(sapply(kest.el, eval,kest.vr2, NULL),
nrow=sqrt(length(kest.el)), byrow=TRUE)
(eigen.analysis(A.mod.b)$lambda1-E$lambda1)/.001
# You could use this strategy to check more rates if you wanted to
# plot the s & e values associated with vital rates
barplot(t(x[,2:3]), beside=TRUE, las=1, xlab="Vital rate",
main=expression(paste("Lesser Kestrels: ", lambda,
"'s sensitivity & elasticity to vital rates")))
abline(h=0)
legend(1,1, rownames(t(x[,2:3])), fill=grey.colors(2))
grid()
## plot effects of changing vital rates from 0 to 1
n<-length(kest.vr)
vr<-seq(0,1,.05) # Sequence of vital rate values
# build storage matrix for lambda obtained at each value of each rate
vrsen<-matrix(NA, nrow=length(vr), ncol=n, dimnames=list(vr, names(kest.vr)))
# for loop works with 1 rate at a time and calculates lambda at each value checked
for (h in 1:n)
{
kest.vr2<-list(b = 0.9321, c0 = 0.3847, ca = 0.925, s0 = 0.3409, sa = 0.7101)
for (i in 1:length(vr))
{
kest.vr2[[h]]<-vr[i]
A<-matrix(sapply(kest.el, eval,kest.vr2 , NULL), nrow=sqrt(length(kest.el)), byrow=TRUE)
vrsen[i,h] <- max(Re(eigen(A)$values))
}
}
# look at the ouptut for the sensitivity analysis
vrsen
# create the plot of lambda (vertical axis) against
# the various values used for each vital rate (x-axis) in a new window
windows()
matplot(rownames(vrsen), vrsen, type='l', lty=1, lwd=2, las=1, col=1:5,
ylab="Lesser kestrel population growth", xlab="Value of vital rate",
main="Effects of changing Lesser Kestrel vital rates")
# plot lambda for original matrix in gray and label that line
abline(h=E$lambda1, col='gray45')
text(0.1,round(E$lambda1,4),bquote(lambda==.(round(E$lambda1,4))),col='gray45')
# plot lambda=1 in black and label that line
abline(h=1)
text(0.1,1.02,bquote(lambda==.(1)))
# build legend
vrn<-expression(b, c0, ca, s0, sa)
legend(.72, .6, vrn, lty=1, lwd=2, col=1:5, cex=1.2)
# put verticals on in corresponding colors for original vital rate values
abline(v=c(.9321,.3847,.925,.3409,.7101),col=1:5,lty=3)
# use the line below if you want to find values for x-y points of interest
# e.g., where lambda=1 for any given vital rate
# the '3' in parentheses dictates that you'll click on 3 x-y coordinates.
# After you do so, R will report the x-y values for those locations.
# I chose 3 because only 'c0', 's0', and 'sa' have lines that cross
# the lambda=1 value
locator(3)