# 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 to make a random draw from, 
# rather than a constant.  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 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)


 # Enter mx data and calculate 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]
	}

# Deterministic projection with Leslie matrix.  This program accomplishes this in two ways
# and compares results as a check.

# load the popbio packages, 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 first year.
Dist_Year[1,] <- N 

# Do the matrix multiplication for 30 years of projection

for(i in 2:31) {
	Dist_Year[i,] <- les.mat %*% Dist_Year[i-1,] 
	}

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

tot <- rep(NA,100) 				# Dummy vector to store 't'-year projection totals.
t <- 31 						# Number of years for the projection.
for(a in 1:100){ 					# Repeat stochastic 't'-year projections a hundred times.
stoc_year <- matrix(NA,nrow=t,ncol=10) 	# Empty matrix of age-class values.
N <- c(97,rep(0,9)) 				# Initial values in each age-class.
stoc_year[1,] <- N 				# Assign 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.	

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.

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.

stoc_year[i,] <- lm %*% stoc_year[(i-1),] # Matrix mult. for next time-step.
} 							# Ends 'i' loop.

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.

# Examine results
tot
ci