#leslie.R
################################################
# DS Wethey 2011-09-12.
# translate leslie.mws to R
################################################
# require (popbio) #needs popbio package
################################################
#R worksheet to accompany
#Brault and Caswell, 1993, Pod-specific demography of killer-whales, Ecology 74: 1444-1454
#and
#Rodrigo Bernal,1998. Demography of the vegetable ivory palm Pytelephas seemanii in Colombia,
#and the impact of seed harvesting. Journal of Applied Ecology 35:64-74
#and other Leslie-Lefkowitch models
#Stage-based models are a modification of the Leslie model. In this case, the assumption is that
# the stage of an organism (size class) is a better indicator of reproductive performance and survival than age is.
# This assumption works well for many plants, and any organism in which size is correlated with survival and fecundity.
#The general form of the model treats the stages as elements of a vector.
#The population projection matrix A (sometimes called a Lefkovich matrix, or, among mathematicians, transition matrix)
# is used to calculate the number of individuals in each stage class at time t+1, based on the numbers at time t.
# n(t+1) = A n(t)
Caswell<- matrix(nrow=4, ncol=4,byrow=T,
c(0,0.0043,0.1132,0,
0.9775,0.9111,0,0,
0,0.0736,0.9534,0,
0,0,0.0452,0.9804 );
Bernal<- matrix(nrow=6,ncol=6,byrow=T,
c(0.7052,0,0,14.8173,14.9770,14.7732,
0.1548,0.6339,0,0,0,0,
0,0.0216,0.9100,0,0,0,
0,0,0.0330,0.9614,0,0,
0,0,0,0.0386,0.9220,0,
0,0,0,0,0.0188,0.9535))
A<-Bernal; A;
n<- nrow(A)
#The elements on the top row are the reproductive outputs.
#The subdiagonal elements represent the probability of transition to the next stage.
#The diagonal elements represent the probability of remaining in a stage.
#here are some questions to think about.
#*In the palm tree model what fraction of seeds remain in the 'seed bank' each year, that is, do not germinate?
#**Of 1000 stage 1 individuals, how many survive to become stage 2?
# If there are currently 1000 stage 2 and 1000 stage 3, how many stage 1 individuals will be present after 1 year?
#***Suppose an individual is in the highest stage. What is the probablility that this individual will still be alive 14 time steps from now?***
# Let's track what happens to a single burst of reproduction, say 1000 seeds or newborns. Each time step is computed by taking the matrix A times the current population vector.
pop<- list(c(1000, 0, 0, 0, 0, 0))
distrib<- list(c( 1, 0, 0, 0, 0, 0));
points<-matrix(ncol=3,nrow=12)
points[1,]<-c(1,sum(pop[[1]]),log10(sum(pop[[1]])))
k<-1
for (i in seq(1,51,5))
{ pop[[i+5]]<-(A%*%A%*%A%*%A%*%A)%*%pop[[i]]
distrib[[i+5]]<-pop[[i+5]]/sum(pop[[i+5]])
points[k+1,]<-c((i+5),sum(pop[[i+5]]),log10(sum(pop[[i+5]])))
k<-k+1}
five_year_rate<- sum(pop[[51]])/sum(pop[[46]]); print(five_year_rate)
annual_rate<- (five_year_rate)^(1/5); print(annual_rate)
par(mfrow=c(1,2))
plot(points[,1],points[,2],type="l", main = "Total population as a function of time" );
plot(points[,1],points[,3],type="l",main = "Log10 of total pop. vs. time");
#Now let's compare this with the eigenvalue - eigenvector analysis.
eigsys<-eigen(A); print(eigsys) #This function returns the eigenvalues and corresponding eigenvectors in order of magnitude
Evec<-eigsys$vectors ;print(Evec)
Eval<-eigsys$values ;print(Eval)
magnitudes<-abs(Eval) ;print(magnitudes)
lambda<-max(magnitudes);print (lambda)
#How does the largest one correspond to the annual rate you calculated earlier?
annual_rate
#If the population grows for many generations with constant transition probabilities given by matrix A,
# the population age or stage vector will approach an eigenvector associated with the dominant eigenvector of A.
# This is the stable age or stage distribution vector.
# Typically the stable stage vector w is scaled so that the sum of the elements is 1.0.
# To scale the stable stage vector, we have to do some fiddling.
# Add up the elements in the vector; then divide all elements in the vector by this sum:
V<-Evec[,1] ;print(V)
#scale this eigenvector so that it sums to 1.0
w<-V/sum(V) ;print(w)
#how does this compare to the last age distribution you calculated?
distrib[[56]]
#Let's check; is A w = (lambda) w?
#Remember this is our definition of an eigenvalue and eigenvector
# - that the eigenvector multiplied by the matrix is the same as the eigenvector multiplied by the eigenvalue.
# If we got the correct eigenvalues and eigenvectors, we should obtain the same answer by multiplying
# the stable age distribution by matrix A or by the eigenvalue:
A %*% w
lambda*w
#sensitivity analysis of Bernal
#the popbio library does all of the eigenvalue/eigenvector calculations and sensitivity/elasticity calculations:
require(popbio)
# if you want the full sensitivity matrix, change to eigen.analysis(A,zero=F)
eigen_analysis<-eigen.analysis(A,zero=T); print(eigen_analysis)
#***Try seeing if you believe that the sensitivities or elasticities are correct, following the pattern illustrated below.***
#We will change one of the entries in the original "A" matrix by 10% and see what percent change there is in lambda.
#Compare this change with the corresponding value in the elasticity or sensitivity matrix.
#Let's first recall the original matrix, eigenvalues, eigenvectors:
old_lambda<-eigen_analysis$lambda1;print(old_lambda)
old_elasticity<-eigen_analysis$elasticities
#Let's change the P2 value by10% (we'll do this by multiplying P2 by 1.1) and see what percent change we get in lambda.
#This is the so-called [2,2] entry of the matrix A, that is, the number that appears in the first row (from top to bottom)
# and second column (from left to right).
#If the author's description of elasticity is correct, then we should get a change in lambda
# that is proportional to the elasticity value for P2.
A[2,2]<-1.1*A[2,2] ;print(A[2,2]) # make the change in the affected entry
print(A)
new_lambda<-eigen.analysis(A)$lambda1; print(new_lambda)
ratio<-new_lambda/old_lambda ;print(ratio)
frac_change<- ratio - 1 ;print(frac_change)
elas<-old_elasticity[2,2] ;print(elas)
#Note that we made a 10 percent change in one element of the matrix (a fractional change of 0.1).
# If you multiply this elasticity value by the fractional change made in the matrix element,
# you should get a number close to the fractional change in lambda.
# If you made a 100% change (a fractional change of 1.0) in the matrix element,
# you should see a fractional change in lambda equal to the elasticity value.
(elas * 0.1) - (frac_change); # close to zero?
#look at eigen analysis before this change in A:
eigen_analysis
#calculate new stable age distribution, reproductive value distribution, etc
eigen.analysis(A)
#The author states on p71 that "an 86% reduction in fecundity of all classes is required to reduce lambda from 1.0589 to 1.0".
# Try making reductions in fecundity and see if your observations confirm this:
#First we will get an original copy of matrix A
A<-Bernal
#Then we define our fractional reduction in fecundity.
fecundity_reduction<-0.86;
#We will now reduce the seed production entries by this amount.
# We will not change A[1,1], because that element of the matrix defines the rate of retention of seeds in the seed bank, not a fecundity.
# Therefore we change only elements 2 to 6 in the top row of the matrix:
for (i in 2:ncol(A))
{A[1,i]<-A[1,i]*(1-fecundity_reduction)}
#rerun the eigen analysis and look at the value of lambda1:
eigen.analysis(A)
#############################################################
#Caswell stable age distribution
rm(list=ls()) # this is equivalent of Maple restart command
Caswell<- matrix(nrow=4, ncol=4,byrow=T,
c(0,0.0043,0.1132,0,
0.9775,0.9111,0,0,
0,0.0736,0.9534,0,
0,0,0.0452,0.9804 ))
A<-Caswell
n<-nrow(A)
#The elements on the top row are the reproductive outputs.
#The subdiagonal elements represent the probability of transition to the next stage.
#The diagonal elements represent the probability of remaining in a stage. here are some questions to think about.
#*In the whale model what tells you that juveniles do not reproduce? that older females do not reproduce?
# Does this remind you of any other creature with a similar life stage pattern?*
#**Of 1000 stage 1 individuals, how many survive to become stage 2?
# If there are currently 1000 stage 2 and 1000 stage 3, how many stage 1 individuals will be present after 1 year?
#***Suppose an individual is in the highest stage.
# What is the probablility that this individual will still be alive 14 time steps from now?***
# Let's track what happens to a single burst of reproduction, say 1000 seeds or newborns.
# Each time step is computed by taking the matrix A times the current population vector.
pop<- list(c(1000, 0, 0, 0))
distrib<- list(c( 1, 0, 0, 0));
points<-matrix(ncol=3,nrow=12)
points[1,]<-c(1,sum(pop[[1]]),log10(sum(pop[[1]])))
k<-1
for (i in seq(1,101,10))
{ pop[[i+10]]<-(A%*%A%*%A%*%A%*%A%*%A%*%A%*%A%*%A%*%A)%*%pop[[i]]
distrib[[i+10]]<-pop[[i+10]]/sum(pop[[i+10]])
points[k+1,]<-c((i+10),sum(pop[[i+10]]),log10(sum(pop[[i+10]])))
k<-k+1}
ten_year_rate<- sum(pop[[101]])/sum(pop[[91]]); print(ten_year_rate)
annual_rate<- (ten_year_rate)^(1/10); print(annual_rate)
par(mfrow=c(1,2))
plot(points[,1],points[,2],type="l", main = "Total population as a function of time" );
plot(points[,1],points[,3],type="l",main = "Log10 of total pop. vs. time");
#Now let's compare this with the eigenvalue - eigenvector analysis.
eigsys<-eigen(A); print(eigsys) #This function returns the eigenvalues and corresponding eigenvectors in order of magnitude
Evec<-eigsys$vectors ;print(Evec)
Eval<-eigsys$values ;print(Eval)
magnitudes<-abs(Eval) ;print(magnitudes)
lambda<-max(magnitudes);print (lambda)
#How does the largest one correspond to the annual rate you calculated earlier?
annual_rate
#If the population grows for many generations with constant transition probabilities given by matrix A,
# the population age or stage vector will approach an eigenvector associated with the dominant eigenvector of A.
# This is the stable age or stage distribution vector.
# Typically the stable stage vector w is scaled so that the sum of the elements is 1.0.
# To scale the stable stage vector, we have to do some fiddling.
# Add up the elements in the vector; then divide all elements in the vector by this sum:
V<-Evec[,1] ;print(V)
#scale this eigenvector so that it sums to 1.0
w<-V/sum(V) ;print(w)
#how does this compare to the last age distribution you calculated?
distrib[[101]]
#sensitivity/elasticity analysis of Caswell
#the popbio library does all of the eigenvalue/eigenvector calculations and sensitivity/elasticity calculations:
require(popbio)
eigen.analysis(A,zero=T)