# BIOL 521: PVA by Stochastic Leslie Matrix Projection
# Multiple matrices approach
# Scott Creel 9/30/2011
# INTRO
# Demographic PVA with Mountain Golden Heather
# from example in Morris & Doak 2002, Quantitative Conservation Biology
# Chapter 6. See their Table 6.7 for data, which come from
# Frost 1990 and Gross et al. 1998.
# Uses a data set for Hudsonia montana, or mountain golden heather. This is a
# plant now restricted to a small portion of N. Carolina (in Pisgah
# National Forest. It is listed as threatened under the ESA, but is considered
# G1 by NatureServe (G = all global populations threatened, 1 = critically
# endangered. The original causes of decline are many (including fire suppression
# and direct losses to human activities. The current range is so restricted
# that the species is at risk of extinction due simple due to small population
# size. Thus a PVA is of interest to establish the magnitude of extinction
# risk. A demographic PVA is of interest so that elasticities (together with
# information on the observed variance in demographic parameteres) can inform
# responses.
# ANNOTATED CODE
# Input the stored data and examine it. It is a series of 4 transition (Leslie) matrices,
# one matrix for each year from 1985 to 1988 each matrix is 6x6 because
# there are 6 STAGE classes. Some elements other than the top row (stage
# specific fertility) and subdiagonal(stage specific probability of transition
# to the next stage) are non-zero in these matrices. These elements represent
# the probability of surviving but staying in the same stage class. In an
# AGE structured matrix, the 'probability of transition to the next stage' is
# the same thing as survival. In a stage structured matrix, survival is
# decomposed into survivors that change stages and survivors that stay in the
# same stage
library(popbio) # load the popbio package
data(hudsonia) # load the Hudsonia data file included with popbio
hudsonia # inspect the data - this is ALWAYS a good idea
### Set up an initial age structure or N0 vector (named 'n' here)
# this N0 vector is realistic for a plant, i.e. highly biased to young
# individuals. Anything that has a large number of offspring that are
# each likely to die -- like seeds -- will have a very flat age pyramid
# i.e. a type III survivorship curve. Recall that this species' life
# history has 'recursive' loops, in which a survivor to remains in the
# same stage class from one year to the next. These are the off-subdiagonal
# elements discussed above. With recurive loops, the stage structure does
# NOT have to show fewer individuals in each successive stage of the life
# history. Individuals stay in the medium stage a long time (probability
# of staying the stage from one year to the next ranges from 0.44 to 0.65.
# So, even if the rate iof individuals enetering the medium stage is low,
# a fair number accumulate.
n<-c(4264, 3,30,16,25,5)
names(n)<-c("seed", "seedlings", "tiny", "small", "medium" , "large")
n
### 1000 replicates of projection for 50 years
# with equal and unequal probabilities for matrix selection
#
# Also, weighting how likely it is that the population will
# experience the demographic conditions captured in each year's
# transition matrix
#
# stoch.projection is a function defined in the popbio package.
# As with any
# function, you can type ?stoch.projection in R to see the arguments for the
# function and its outputs. You can do the same in RStudio by
# putting the cursor within the function and hitting the F1 key.
# In this example, notice in the code below that there
# is no argument for the number of years projected. 50 is the default, so if you
# omit "tmax = [any number]", you get 50 years.
reps = 1000
x.eq<-stoch.projection(hudsonia, n, nreps=reps) #each matrix has equal probability of occuring
x.uneq<-stoch.projection(hudsonia, n, nreps=reps, prob=c(.2,.2,.2,.4)) #last matrix more likely than the others
# plot frequency distribution of the population size in year 50 for the 1000
# replicates of each simulation, with a dashed vertical line indicating
# N0, initial population size. Because this is a stochastic projection, you
# won't get identical results each time. If we increase the number of replicates
# (5000 is the maximum allowed by the stoch.projection function), the results
# will be more stable. Using a small number of replicates and comparing the differences
# among runs is actually an indirect method of examining ENVIRONMENTAL STOCHASTICITY.
par(mfrow =c(2,1)) # set up the plot area to display two plots, one above the other (2 rows, 1 column)
hist(apply(x.eq, 1, sum), xlim=c(0,5000), ylim=c(0,200), col="green",
breaks=seq(0,5000, 100), xlab="Final population size at t=50",
main='Projection of stochastic growth for Hudsonia
using equal and unequal probabilities of drawing each matrix')
## add initial pop size to the plot
sum(n)
abline(v=sum(n), lty=3)
## use transparency for overlapping distributions - may not work on all systems
hist(apply(x.uneq, 1, sum), xlim=c(0,5000), ylim=c(0,200), col = rgb(0, 0, 1, 0.2),
breaks=seq(0,10000, 100),xlab="Final population size at t=50", main = '' )
## initial pop size to the plot
sum(n)
abline(v=sum(n), lty=3)
### Extensions - once you have run the model, do the following:
# 1. Calculate the probability of pseudo-exinction, defined as dropping
# below ext.threshold, initially set to 500 individuals in this example.
N.final = apply(x.eq, 1,sum) # sum the age classes to give N at final time step for each rep
ext.threshold = 500 # set the definition of pseudoextinction
ext = 0 # set the counter for number of extinctions to zero before tallying extinctions
# loop to tally up the number of cases where N.final is < ext.threshold
for (i in 1:length(N.final))
{if (N.final[i] < ext.threshold) ext = ext+1
}
ext # number of cases where N.final is < ext.threshold
prob.ext = ext/reps
prob.ext # proportion of cases where N.final is < ext.threshold, i.e. Prob{pseudoextinction}
# 2. Repeat step one for several different extinction thresholds. Is the relationship
# between prob.ext and ext.threshold linear?
N.final = apply(x.eq, 1,sum) #could use x.eq or x.uneq
ext.threshold = 0
increment = 100
ext=matrix(0,10,1)
prob.ext=matrix(0,10,1)
for (i in 1:10)
{ext.threshold = (ext.threshold + increment)
ext[i] = 0
for (j in 1:length(N.final))
{if (N.final[j] < ext.threshold) ext = ext+1
}
ext[i]
prob.ext[i] = ext[i]/reps
prob.ext[i]
}
ext
prob.ext
thresholds = matrix(seq(0,900, by = 100),10,1)
par(mfrow = c(1,1)) #reset the plot area to show one plot
plot(thresholds, prob.ext, xlab = 'Pseudo-extinction threshold',
ylab = 'Probability of Extinction')
# fit linear model
lin.mod=lm(prob.ext~thresholds)
# get coefficient estimates (intercept and slope)from linear model
out=matrix(NA,1,2)
out[1:2]=coef(lin.mod)[1:2]
# add fitted model to plot
abline(a = out[,1], b = out[,2], lty = 1)
# 2. Incorporate density dependence and evaluate the effect on population growth and
# pseudo-extinction. This simulates the effect of systematic (not stochastic) changes
# in environmental conditions.
# ?stoch.projection will let you examine all of the arguments of this function.
# Nmax is one of the arguments, and it allows you to add density dependence
# to the projection. This is 'ceiling'density dependence, the simplest possible version
# of DD growth, where the population grows exponentially until it hits the ceiling
# defined by Nmax, and then abruptly drops to r <= 0 or lambda <= 1.
# 3. Examine the effect of different initial age structures and population sizes
# on the probability of pseudo-extinction. This lets you consider how the
# demographic circumstances affect demographic stochasticity and thus extinction risk.
#