# BIOL 521: PVA by Stochastic Leslie Matrix Projection
# Approach using matrix entries that specify a
# mean and variance (and a specific distribution)
# for each demographic variable in the matrix
# Scott Creel 11/17/2011
# First, clear memory and load the popbio package (install it firstif you have not already done so)
rm(list=ls())
library(popbio)
# INTRODUCTION
# This exercise builds a stochastic Leslie matrix to estimate extinction risk for
# African wild dogs, using the data from Tanzania's Selous Game Reserve that we used
# in class to build a life table and estimate lambda.
#
# In the previous exercise you saw the 'multiple matrices' approach to
# stochastic projection... This approach assembles a set of Leslie matrices and randomly
# draws one entire matrix for each step of projection. This approach assumes that all of the
# demographic parameters covary -- they ALL change at once.
# This assumption might not be correct, so you might want to do the projection in a manner that allows the
# demographic parameters to vary independently of one another. To do this,
# you must make random draws from a distribution for each demographic parameter,
# rather than randomly drawing whole matrices.
#
# All that this requires is:
# 1. Create a Leslie matrix in which each of the demographic parameters is specified by
# a distribution rather than a constant. You do this by specifying a distribution, a mean and a variance.
# There are other logical possibilities, but we'll use a binomial distribution for survival
# rates and a Poisson distribution for fecundities.
# 2. Making a random draw from those distributions at each time step to get a 'stochastic Leslie matrix'.
# 3. Projecting next year's population (tracked as a vector with the number
# of individuals in each age class) with the stochastic Leslie matrix, just the same way you
# would with a simple Leslie matrix with the mean for each demographic variable.
# PART 1: BASIC DEMOGRAPHIC CALCULATIONS, SURVIVORSHIP CURVE.
# These are the African wild dog data we used to build a life table in class.
# Enter lx data (survivorship from birth) and examine survivorship curve
lx <- c(1,0.75,0.63,0.49,0.37,0.29,0.13,0.09,0.02,0.02,0.02)
log.lx <- log(lx)
# create age variable x that corresponds to the lx values above
x <- 0:10
plot(x,log.lx,main="Survivorship Curve",xlab="Age Class - x",
ylab="Survivorship - log(lx)",pch=16)
# add a line for perfect type II survivorship, just for visual comparison
# type II survivorship is the case with a constant survival rate (sx)
# for all age classes, which yields a straight line
# on a survivorship curve with a log-Y axis
lines(x,(-3.9/10*x),type="l",lty=2)
# Enter mx data (fecundity) and calculate approximate growth rate
# These values differ from the ones in the life table in class, because they
# are for a post-birth pulse annual census approach
mx <- c(0,0,0.1,0.28,0.59,1.12,3.99,2.14,3.32,2.66,0)
contrib.growth <- lx*mx # lx*mx for each age class
Ro <- sum(contrib.growth) # summed to get Ro
x.lxmx <- x*contrib.growth # x*lx*mx for each age class
T <- sum(x.lxmx)/Ro # summed and divided by Ro to get generation time
r.app <- log(Ro)/T # and then used to get approximate value of r
# Look at the results
Ro
T
r.app
# Use euler equation to check if the approximate solution for r is reasonable.
# If the sum is close to one, the approximate solution is fairly accurate.
# You Can get exact solution by altering the estimated value of r until the sum
# equals one, but can much more easily get exact lambda [see below] as the dominant
# Eigenvalue of the Leslie matrix (and then r = ln(lambda))
sum(exp(-r.app*x)*lx*mx) # pretty close to one, approximate solution is not terrible for these data
# PART 2: DETERMINISTIC LESLIE PROJECTION
# convert lx data to sx values - you may have to think a bit to see why this works
# (making the lx and sx columns on a piece of paper helps a lot)
Sx <- c(73/97,61/73,48/61,36/48,28/36,13/28,9/13,2/9,2/2,2/2,0/2)
#convert mx to Fx (here for a post-birth pulse count)
Fx <- mx*Sx[1] # fecundity values
S <- Sx[-1] # dropping the first value
F <- Fx[-1]
# create an empty matrix (with all zeros) of the correct dimensions.
les.mat <- matrix(rep(0,10*10),nrow=10)
# Fill in relevant values for the non-zero entries
# -- row 1 is fecundities
les.mat[1,] <- F
# -- and the subdiagonal is annual survival rates, using a loop to efficiently assign them to
# the correct positions within the matrix
for(i in 1:9){
les.mat[(i+1),i] <- S[i]
}
#examine the completed projection matrix: Fx values in the top row, Sx in subdiagonal, zeros elsewhere
les.mat
# Deterministic projection with Leslie matrix. This script does this in two ways
# and compares results as a check. First you'll do it for yourself, and then use
# the popbio package, which has functions that make it much easier
# Input a vector with the initial distribution of individuals
# in across age classes. In this example, 97 newborns only.
N <- c(97,0,0,0,0,0,0,0,0,0)
# Create an empty matrix [with NA values] to store the set of projected population vectors
# (as 10 columns) for 30 years of projection into the future (each as row), starting with population
# defined by N above.
Dist_Year <- matrix(NA,nrow=31,ncol=10)
# Set the initial population in the first year.
Dist_Year[1,] <- N
# Do the matrix multiplication for 30 years of projection, years 2 to 31
# matrix multiplication uses %*% rather than *
for(i in 2:31) {
Dist_Year[i,] <- les.mat %*% Dist_Year[i-1,]
}
# Calculate and examine total size at each time step. In the apply() function, the sum argument
# species to calculate a sum, the 1 specifes that rows should be summed (2 would specify to sum
# columns), and the Dist_Year indicates what object to do all this on. This is a very old-school
# line of code, and there are MANY other ways to do this, either using the sum() command or using the
# summarise command in dplyr. This code is using base R only.
proj <- apply(Dist_Year,1,sum)
proj
# Double-check using the Popbio package function 'pop.projection'.
p <- pop.projection(les.mat,N,31)
p
# Examine the results. Note that results for popualtion size agree exactly, and that the
# pop.projection function also provides a lot of other useful output,
# including lambda (dominant eigenvalue of the matrix) and the stable
# age distribution (R eigenvector of the matrix). The stable age distribution
# output is output as the proportion of the population in each age class.
# Popbio includes a function eigen.analysis that will calculate the above information, plus
# sensitivities and elasticities.
# Explore. Remember that you can press F1 with the cursor in a function to get help on it.
eigen.analysis(les.mat)
# PART 3: STOCHASTIC LESLIE PROJECTION
# This uses 4 nested loops. The j and k loops are innermost and draw stochastic fecundity and survival.
# The t loop is next, and steps across 30 years of projection from an initial N vector
# the a loop is outermost and replicates the 30 years of projection for 100 iterations.
tot <- rep(NA,100) # Dummy vector to store total population size for 100 replicates
t <- 31 # Number of years for the projection.
for(a in 1:100){ # Repeat stochastic 't'-year projections 100 times.
stoc_year <- matrix(NA,nrow=t,ncol=10) # Dummy matrix of age-class values. One row per time step. Columns are # individuals in each age class.
N <- c(97,rep(0,9)) # Initial values in each age-class.
stoc_year[1,] <- N # Assign these values to first row in age-class matrix.
for(i in 2:t){ # Begin loop for each projection over t yeats
lm <- les.mat # Set up a temporary Leslie matrix for each iteration
# with the correct mean fecundities and survival rates
for(j in 1:10){ # Randomly draw fecundities for each class.
lm[1,j] <- rpois(1,les.mat[1,j]) # Mean of Poisson is fecundity value.
} # Ends the 'j' loop that draws fecundities.
for(k in 1:9){ # Randomly draw survival probabilities for each class.
n <- stoc_year[(i-1),k] # n is number of individuals to calculate live/die for.
lm[(k+1),k] <- ifelse(n>1,rbinom(1,size=round(n),p=les.mat[(k+1),k])/n,0)
# need ifelse statement to deal with the possibility
# that there are no individuals in that age class.
} # Ends the 'k' loop to draw survival rates.
stoc_year[i,] <- lm %*% stoc_year[(i-1),] # Matrix multiplication for next time-step.
} # Ends the 'i' loop projecting across time.
tot[a] <- sum(stoc_year[t,])
} # Ends the 'a' loop of one hundred iterations of the projections
ci <- quantile(tot,probs=c(0.025,0.975)) # Find 95% confidence intervals on final population size.
# Examine final population size for each of 100 iterations, and the 95% confidence interval on final population size
# This is stochastic so I can't tell you exactly what numbers you'll get, but your confidence interval
# on population size 30 years from now probably does not overlap zero.
tot
ci