# Matrix analysis of thistle population adapted from appendix in following article
# Ezard et al. 2010. Matrix models for a changeable world:
# the importance of transient dynamics in population management.
# J. Applied Ecology 47:515-523
require(popbio)
## Stages
# 1. SB = seeds in seed bank
# 2. S = small rosettes
# 3. M = medium rosettes
# 4. L = large rosettes
## Vital-rate definitions
# sgma1 Survival of seed in SB
# sgma2 Survival of S (small rosettes)
# sgma3 Survival of M (medium rosettes)
# sgma4 Survival of L (large rosettes)
# gam3 Growth of establishing seed to M
# gam32 Growth of surviving, not-bolting S to M
# gam42 Growth of surviving, not-bolting S to L
# gam43 Growth of surviving, not-bolting M to L
# beta2 Bolting of surviving S
# beta3 Bolting of surviving M
# beta4 Bolting of surviving L
# pi2 Potential seed production by S
# pi3 Potential seed production by M
# pi4 Potential seed production by L
# phi % seeds escaping from floral herbivory
# nu New seed entering SB
# eps New seed establishing seedling
# eps1 Seed from SB establishing seedling
# vital rate values stored as list
thistle.vr <- list(sgma1=0.2597, sgma2=0.2182, sgma3=0.5633, sgma4=0.7474,
gam3=0.1627, gam32=0.8028, gam42=0.1268, gam43=0.3824,
beta2=0.0180, beta3=0.3596, beta4=1,
pi2=100, pi3=195, pi4=1473,
phi=1, nu=0.2333, eps=0.03, eps1=0.03)
# matrix elements stored as expressions by row of matrix
thistle.el <- expression(
# row 1 of matrix
sgma1, sgma2*beta2*pi2*phi*nu, sgma3*beta3*pi3*phi*nu, sgma4*beta4*pi4*phi*nu,
# row 2 of matrix
eps1*(1-gam3), sgma2*(1-beta2)*(1-gam32-gam42) + sgma2*beta2*pi2*phi*eps*(1-gam3),
sgma3*beta3*pi3*phi*eps*(1-gam3), sgma4*beta4*pi4*phi*eps*(1-gam3),
# row 3 of matrix
eps1*gam3, sgma2*(1-beta2)*gam32 + sgma2*beta2*pi2*phi*eps*gam3,
sgma3*(1-beta3)*(1-gam43) + sgma3*beta3*pi3*phi*eps*gam3, sgma4*beta4*pi4*phi*eps*gam3,
# row 4 of matrix
0, sgma2*(1-beta2)*gam42, sgma3*(1-beta3)*gam43, sgma4*(1-beta4))
# evaluate matrix element values by applying vital rate values
A=matrix(sapply(thistle.el,eval,thistle.vr,NULL),nrow=sqrt(length(thistle.el)), byrow=TRUE)
rownames(A)=c("seeds","sm.ros","med.ros","lg.ros")
colnames(A)=rownames(A)
# print out matrix
A
# Conduct eigen analysis to obtain various asymptotic quantities of interest
E=eigen.analysis(A)
E
# get the s & e values associated with vital rates
VitSens<-vitalsens(thistle.el, thistle.vr)
VitSens
# plot the elasticity values (3rd column)
bp=barplot(t(VitSens[,3]), las=1, xlab="Vital rates",
ylab=expression(paste("Elasticity of ", lambda ,
" to vital rates")),
names.arg=c("sgma1","sgma2","sgma3","sgma4","gam3",
"gam32","gam42","gam43","beta2","beta3",
"beta4","pi2", "pi3","pi4", "phi",
"nu", "eps", "eps1"), las = 3)
abline(h=0)