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