# 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 and Schulenberg Prairies. # 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 (B) and Schulenberg (S) Prairies divided into # between ecotype and between population crosses (bebp), # within ecotype between population (webp) and # within ecotype within population (wewp) # 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) # this version uses fecundity calculated from seed germination rates from cross records # (M. Bowles, unpublished data) library(popbio) stages <- c("sdlg", "juv1", "juv2", "juv3", "flow") # ------------------------------------------------------------------------------- # This section calculates lambda and variance for each matrix # The matrix names are below indicating their population (B or S) and crosses (bebp, webp, wewp) # A1 = Bbebp # A2 = Bwebp # A3 = Bwewp # A4 = Sbebp # A5 = Swebp # A6 = Swewp # Transition matrices as vectors A1 <- matrix(c(0, 0, 0, 0, 9.7416, 0.95, 0.896648045, 0, 0, 0, 0, 0.01396648, 0.777777778, 0, 0, 0, 0.005586592, 0.148148148, 0.75, 0, 0, 0, 0, 0.25, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) A2 <- matrix(c(0, 0, 0, 0, 15.4011, 0.9, 0.873873874, 0.5, 0, 0, 0, 0.009009009, 0.49, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0.01, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) A3 <- matrix(c(0, 0, 0, 0, 5.9268, 0.794117647, 0.658333333, 0.093023256, 0.035087719, 0.125, 0.029411765, 0.195833333, 0.720930233, 0.245614035, 0.25, 0, 0.008333333, 0.093023256, 0.631578947, 0.375, 0, 0.0125, 0.023255814, 0.052631579, 0.25), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) A4 <- matrix(c(0, 0, 0, 0, 9.4083, 0.914893617, 0.784090909, 0.166666667, 0, 0, 0, 0.015151515, 0.5, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0.01, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) A5 <- matrix(c(0, 0, 0, 0, 20.3511, 0.96, 0.747368421, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0.01, 0, 0, 0, 0, 0, 0.01, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) A6 <- matrix(c(0, 0, 0, 0, 8.25, 0.857142857, 0.74501992, 0.112094395, 0.012048193, 0.071428571, 0, 0.087649402, 0.675516224, 0.21686747, 0.464285714, 0, 0.011952191, 0.061946903, 0.530120482, 0.178571429, 0, 0.003984064, 0.041297935, 0.072289157, 0.214285714), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) # Transition matrices A1 A2 A3 A4 A5 A6 # calculation of sensitivity sensb<-sensitivity(A1, zero=FALSE) sensc<-sensitivity(A2, zero=FALSE) sensa<-sensitivity(A3, zero=FALSE) sensd<-sensitivity(A4, zero=FALSE) sense<-sensitivity(A5, zero=FALSE) sensf<-sensitivity(A6, zero=FALSE) # conversion of sensitivity matrix to a vector sensbvec <- c(sensb) senscvec <- c(sensc) sensavec <- c(sensa) sensdvec <- c(sensd) sensevec <- c(sense) sensfvec <- c(sensf) # Sensitivity matrices sensb sensc sensa sensd sense sensf # Standard error matrices sea <- matrix(c(0, 0, 0, 0, 0.71706, 0.13788775, 0.030477569, 0.002657349, 0.000523706, 0, 0.020427815, 0.023218093, 0.021122522, 0.003665945, 0, 0, 0.000294605, 0.002084196, 0.005290504, 0, 0, 0.000445145, 0.000565394, 0.000696868, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) seb <- matrix(c(0, 0, 0, 0, 1.17858, 0.222222222, 0.044078069, 0, 0, 0, 0, 0.006086388, 0.027217537, 0, 0, 0, 0.000230971, 0.004241694, 0, 0, 0, 0, 0, 0, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) sec <- matrix(c(0, 0, 0, 0, 1.863296, 0.0526315789473684, 0.094599736, 0, 0, 0, 0, 0.006826785, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) sed <- matrix(c(0, 0, 0, 0, 1.13826, 0.03730628, 0.068517608, 0, 0, 0, 0, 0.004901166, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) see <- matrix(c(0, 0, 0, 0, 2.46216, 0.166666667, 0.087230966, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 5, byrow = TRUE, dimnames = list(stages, stages)) sef <- matrix(c(0, 0, 0, 0, 0.99812, 0.076923077, 0.062044774, 0.009729863, 0.00119356, 0, 0, 0.024293266, 0.058830027, 0.018078069, 0, 0, 0.000723597, 0.003951005, 0.024650671, 0, 0, 0.000243694, 0.002775834, 0.004797416, 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... cova<-diag(unlist(c(sea))^2) covb<-diag(unlist(c(seb))^2) covc<-diag(unlist(c(sec))^2) covd<-diag(unlist(c(sed))^2) cove<-diag(unlist(c(see))^2) covf<-diag(unlist(c(sef))^2) ## calculate the variance of lambda from step 8 in Slakski 2007 vara<-t(sensavec) %*% ( cova %*% sensavec) varb<-t(sensbvec) %*% ( covb %*% sensbvec) varc<-t(senscvec) %*% ( covc %*% senscvec) vard<-t(sensdvec) %*% ( covd %*% sensdvec) vare<-t(sensevec) %*% ( cove %*% sensevec) varf<-t(sensfvec) %*% ( covf %*% sensfvec) # ------------------------------------------------------------------ # this section prints the output (transition matrix, se matrix, sensitivity matrix, lambda, variance, and 95%CI) # the matrices print("Biesecker between ecotype between population") A1 seb sensb print("Biesecker within ecotype between population") A2 sec sensc print("Biesecker within ecotype within population") A3 sea sensa print("Schulenberg between ecotype between population") A4 sed sensd print("Schulenberg within ecotype between population") A5 see sense print("Schulenberg within ecotype within population") A6 sef sensf # lambda and CI print("Biesecker between ecotype between population") 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("Biesecker within ecotype between population") lambda(A2) varc # +/- CI 1.96*sqrt(varc) # and the confidence intervals lambda(A2) - 1.96*sqrt(varc) lambda(A2) + 1.96*sqrt(varc) print("Biesecker within ecotype within population") lambda(A3) vara # +/- CI 1.96*sqrt(vara) # and the confidence intervals lambda(A3) - 1.96*sqrt(vara) lambda(A3) + 1.96*sqrt(vara) print("Schulenberg between ecotype between population") lambda(A4) vard # +/- CI 1.96*sqrt(vard) # and the confidence intervals lambda(A4) - 1.96*sqrt(vard) lambda(A4) + 1.96*sqrt(vard) print("Schulenberg within ecotype between population") lambda(A5) vare # +/- CI 1.96*sqrt(vare) # and the confidence intervals lambda(A5) - 1.96*sqrt(vare) lambda(A5) + 1.96*sqrt(vare) print("Schulenberg within ecotype within population") lambda(A6) varf # +/- CI 1.96*sqrt(varf) # and the confidence intervals lambda(A6) - 1.96*sqrt(varf) lambda(A6) + 1.96*sqrt(varf)