# 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 parameter
# Scott Creel 11/17/2011
# INTRO
#
# In the previous exercise you saw the 'multiple matrices' approach to
# stochastic projection... assemble a set of Leslie matrices and randomly
# draw one for each step of projection. This approach assumes that all of the
# demographic parameters covary, i.e., the all change at once. This might not
# be true, 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 creating a Leslie matrix in which each of the
# demographic parameters is specified by a distribution rather than a constant, and making a random draw
# from that distribution at each time step. You do this by specifying a distribution, a mean and a variance.
# There are more complex possibilities, but we'll use a binomial distribution for survival
# rates and a normal distribution (truncated to be non-negative) for fecundities.
# PART 1: BASIC DEMOGRAPHIC CALCULATIONS, SURVIVORSHIP CURVE
# Enter lx data 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)
x <- 0:10
plot(x,log.lx,main="Survivorship Curve",xlab="Age Class - x",
ylab="Survivorship - log(lx)",pch=16)
lines(x,(-3.9/10*x),type="l",lty=2) # add a kine for perfect type II survivorship, just for visual comparison
# Enter mx data and calculate approximate growth rate
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
Ro <- sum(contrib.growth)
x.lxmx <- x*contrib.growth
T <- sum(x.lxmx)/Ro
r.app <- log(Ro)/T
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.
sum(exp(-r.app*x)*lx*mx)
# PART 2: DETERMINISTIC LESLIE PROJECTION
# convert lx data to sx values
Sx <- c(73/97,61/73,48/61,36/48,28/36,13/28,9/13,2/9,2/2,2/2,0/2) # Survival
#convert mx to Fx (here for a post-birth pulse count)
Fx <- mx*Sx[1] # Fecundity values
S <- Sx[-1]
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 assign them to
# the correct positions within the matrix
for(i in 1:9){
les.mat[(i+1),i] <- S[i]
}
les.mat #examine the completed projection matrix
# Deterministic projection with Leslie matrix. This program accomplishes this in two ways
# and compares results as a check.
# load the popbio package, which has functions that do all of this directly
library(popbio)
# 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 to store the set of projected population vectors
# for 30 years of projection into the future, 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
for(i in 2:31) {
Dist_Year[i,] <- les.mat %*% Dist_Year[i-1,] # matrix multiplication uses %*%
}
# Calculate and examine total size at each time step
proj <- apply(Dist_Year,1,sum)
proj
# Double-check using the Popbio package function 'pop.projection'.
p <- pop.projection(les.mat,N,31)
# 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).
p
# Note also that Popbio includes a function eigen.analysis that will
# give you the elasticities. Explore.
# PART 3: STOCHASTIC LESLIE PROJECTION
# This uses 4 nested loops. The i and j 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.
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 projections.
lm <- les.mat # Setup 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 true fecun. value.
} # Ends the 'j' loop to draw fecundities.
for(k in 1:9){ # Randomly draw survival probabilities for each class.
n <- stoc_year[(i-1),k]
lm[(k+1),k] <- ifelse(n>1,rbinom(1,size=round(n),p=les.mat[(k+1),k])/n,0)
} # Ends the 'k' loop to draw survival rates.
stoc_year[i,] <- lm %*% stoc_year[(i-1),] # Matrix mult. for next time-step.
} # Ends the 'i' loop projecting across time.
tot[a] <- sum(stoc_year[t,])
} # Ends the 'a' loop; hundred individual projections are done.
ci <- quantile(tot,probs=c(0.03,0.97)) # Find 94% confidence intervals on pop size.
# Examine results
tot
ci