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