# 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