# 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)