# R code for calculating lambda and 95% confidence interval for demographic transition matrices

# R version 2.10.0 (2009-10-26)
# Copyright (C) 2009 The R Foundation for Statistical Computing
# ISBN 3-900051-07-0

# This R code was used to calculate population growth rate (lambda) and 95% confidence limits
# for transition matrices from plants established in restorations at Biesecker Prairie.
# See the following publication for details:
# Bowles, M.L., J. L. McBride, and T.J. Bell.  In Review.  Long-term restoration processes affect 
# population growth of the Federal Threatened Mead's Milkweed (Asclepias meadii). Ecosphere.

# Mead's milkweed demographic data for 1994-2011 from 
# Biesecker Prairie divided into burn and unburn plots
# code revised from Hawaiian hawk example in Skalski 2007 (Skalski, J. R., 
# J. J. Millspaugh, P. Dillingham, R. A. Buchanan. 2007. 
# Calculating the variance of the finite rate of population change from a 
# matrix model in Mathematica. Environmental Modelling & Software, 22: 359-364)

library(popbio) 

stages <- c("sdlg", "juv1", "juv2", "juv3", "flow")

# burn matrix, 0.0000, 0.0000, 0.0000, 0.0000, 18.9750, 0.9028, 0.7652, 0.0864, 0.1111, 0.0000, 0.0000, 0.0652, 0.7284, 0.3889, 0.3333, 0.0000, 0.0130, 0.0988, 0.3889, 0.3333, 0.0000, 0.0043, 0.0247, 0.1111, 0.3333
# unburn matrix, 0.0000, 0.0000, 0.0000, 0.0000, 11.9167, 0.9726, 0.8117, 0.1506, 0.0000, 0.0000, 0.0274, 0.0568, 0.6807, 0.2286, 0.0000, 0.0000, 0.0000, 0.0723, 0.6857, 0.5000, 0.0000, 0.0000, 0.0181, 0.0571, 0.2500
# burn se, 0, 0, 0, 0, 8.116466673, 0, 0.08509847, 0.003520082, 0, 0, 0, 0.00725271, 0.029669265, 0, 0, 0, 0.001450542, 0.004022951, 0, 0, 0, 0.000483514, 0.001005738, 0, 0
# unburn se, 0, 0, 0, 0, 4.161926316, 0, 0.037855718, 0.005067855, 0, 0, 0, 0.002650729, 0.022906704, 0.003565062, 0, 0, 0, 0.00243257, 0.010695187, 0, 0, 0, 0.000608143, 0.000891266, 0


# -------------------------------------------------------------------------------
# This section calculates lambda and variance for each matrix
# The matrix names are below indicating their treatment
# A1 = burn
# A2 = unburn

# Transition matrices as vectors
A1 <- matrix(c(0.0000, 0.0000, 0.0000, 0.0000, 18.9750, 0.9028, 0.7652, 0.0864, 0.1111, 0.0000, 0.0000, 0.0652, 0.7284, 0.3889, 0.3333, 0.0000, 0.0130, 0.0988, 0.3889, 0.3333, 0.0000, 0.0043, 0.0247, 0.1111, 0.3333), nrow = 5, byrow = TRUE, dimnames = list(stages, stages))
A2 <- matrix(c(0.0000, 0.0000, 0.0000, 0.0000, 11.9167, 0.9726, 0.8117, 0.1506, 0.0000, 0.0000, 0.0274, 0.0568, 0.6807, 0.2286, 0.0000, 0.0000, 0.0000, 0.0723, 0.6857, 0.5000, 0.0000, 0.0000, 0.0181, 0.0571, 0.2500), nrow = 5, byrow = TRUE, dimnames = list(stages, stages))

# Transition matrices
A1
A2

# calculation of sensitivity 
sensb<-sensitivity(A1, zero=FALSE)
sensu<-sensitivity(A2, zero=FALSE)

# conversion of sensitivity matrix to a vector
sensbvec <- c(sensb)
sensuvec <- c(sensu)

# Sensitivity matrices
sensb
sensu

# Standard error matrices
seb <- matrix(c(0, 0, 0, 0, 8.116466673, 0, 0.08509847, 0.003520082, 0, 0, 0, 0.00725271, 0.029669265, 0, 0, 0, 0.001450542, 0.004022951, 0, 0, 0, 0.000483514, 0.001005738, 0, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages))
seu <- matrix(c(0, 0, 0, 0, 4.161926316, 0, 0.037855718, 0.005067855, 0, 0, 0, 0.002650729, 0.022906704, 0.003565062, 0, 0, 0, 0.00243257, 0.010695187, 0, 0, 0, 0.000608143, 0.000891266, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages))

# Covariance matrix.  In Step 7 in Slakski et al 2007 ("Calculating the variance of the finite rate of population change from a matrix model in Mathematica") they use the square of the standard errors, so we do the same... 
covb<-diag(unlist(c(seb))^2)
covu<-diag(unlist(c(seu))^2)

## calculate the variance of lambda from step 8 in Slakski 2007
varb<-t(sensbvec)  %*%  ( covb %*% sensbvec)
varu<-t(sensuvec)  %*%  ( covu %*% sensuvec)

# ------------------------------------------------------------------
# this section prints the output (transition matrix, se matrix, sensitivity matrix, lambda, variance, and 95%CI)

# the matrices
print("burn")
A1
seb
sensb
print("unburn")
A2
seu
sensu

# lambda and CI
print("burn")
lambda(A1) 
varb
# +/- CI
1.96*sqrt(varb)
# and the confidence intervals 
lambda(A1) - 1.96*sqrt(varb) 
lambda(A1) + 1.96*sqrt(varb) 
print("unburn")
lambda(A2) 
varu
# +/- CI
1.96*sqrt(varu)
# and the confidence intervals 
lambda(A2) - 1.96*sqrt(varu) 
lambda(A2) + 1.96*sqrt(varu)