This week at the Davis R Users’ Group we had a great presentation by Kara Moore O’Leary on using the popbio package to examine rare plant demographics. The following is her script run through knitr. You can download the original script and associated data here. Find out more about Kara and her work at her website here

Demographic analysis using the popbio library and some other fun stuff

A population viability type analysis for a rare herbacious perennial plant, Penstemon albomarginatus, for its only remaining California population

Objectives:

1. Estimating the population growth rate – deterministic and stochastic methods
2. Sensitivity and elasticity analysis – which transitions/rate are most variable, sensitive to change
3. Projecting the population stochastically for different scenarios, observing variation in growth rate
4. Quasi extinction probability
5. Estimating vital rates from real data, which will never be as good/much the data that you would like to have for a PVA – some therapy for me, and maybe you too

Load popbio and some other helpful libraries. I use the dropdown menu in RStudio to install libraries, it seems to work more consistently for my mac than scripting the installs, sorry.
Am I missing something basic about installing packages?

library(popbio)
library(plyr)
library(reshape)

I have already done a bit of magic on my data get here, most or all of which is suggested in the Morris and Doak 2005 book Quantitative population biology. That book was an excellent guide to this kind of analysis. If you want to go deep (and deeper) into this kind of analysis you should get borrow or steal Morris and Doak 2005 and Caswell 2002 from the start! They are both really excellent and straightforward and you’ll be able to answer a lot of the questions that you have along the way. Most of the functions in popbio are based on functions in these two books. Many of the popbio examples are on their data and generate the same graphs etc.

What you need for a class/stage structured demographic model:

• a bunch of individuals (maybe in different populations)
• annual survival rates
• annual class or stage transition rates
• annual fecundity, ie. probability of contribution to the juvenile class in the next year
• these data for a lot of years, preferably 10 or more (but we do what we can with less)

The crux for plants is that its challenging to 1) count all seeds produced annually and 2) know how many seeds really yeilds a juvenile in any year. I am doing a bit of work on these issues from limited field data. If you have a well behaved penguin that you can collar and track and you know produces 2 live juveniles each season, your PVA might be a bit simpler.

Our field data had surivival for each year, mean plant diameter, and inflorescence count classes. From this I found the median inflorescence number for each class MEDIAN_INLF I assigned each plant to CLASS based on its xDIAM_cm. I convinced myself that this had biological meaning by looking at the relationship between size class and survival. Whatever classes or stages you use, you should be confident that they are meaningful for your study species. OR use an integral projection model (IPM) instead. These allow you to use continuous variables like size or age rather than classes. nice abbreviated dataset

andre <- read.csv("karadat.csv")
str(andre)
## 'data.frame':    1221 obs. of  4 variables:
##  $ PLANT_UNQ: Factor w/ 395 levels "1_703","1_704",..: 238 240 249 253 258 263 250 255 259 239 ...
##  $ YEAR     : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS    : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ SEEDS    : num  71.8 0 71.8 0 71.8 ...
head2(andre)
##      PLANT_UNQ YEAR CLASS SEEDS
## 1          A_1 1994    A1 71.75
## 2          A_3 1994    A1  0.00
## 3         B_13 1994    A1 71.75
## .            .    .     .     .
## 1221    6_772A 2012     J  0.00

Let’s look at the stages/classes

levels(andre$CLASS)
## [1] "A1"   "A2"   "A3"   "A4"   "dead" "J"
# dead. After that it can be omitted.

1. ESTIMATING FECUNDITY

The big challenge for plants = estimating seeds/indiv –> juveniles produced the next year based on a few fruit and seed counts and a lot of inflorescence class data

Part 1: How to estimate seed production per plant?

Multiply seeds/infl from field observations by inflorescence number. Ideally you would have seed counts for each plant, but in the absence of those nearly impossible data, I’m using a few seeds/fruit counts * fruits/infl counts from 2011 and 2012.

Part 2: How many seeds makes an emergent juvenile?

I’ll add some notes on how I made these estimates at the end of this script, but let’s start today with a dataset ready to go for popbio I’m sure that every PVA has a hurdle or five to leap before you have the vital rate data ready to go. I’ve added estimated seeds produced per individual per year as a column in this simple dataset more on how I generated that later if we are interested

str(andre)
## 'data.frame':    1221 obs. of  4 variables:
##  $ PLANT_UNQ: Factor w/ 395 levels "1_703","1_704",..: 238 240 249 253 258 263 250 255 259 239 ...
##  $ YEAR     : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS    : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ SEEDS    : num  71.8 0 71.8 0 71.8 ...

2. PICK A STARTING POPULATION VECTOR

this is the # of individuals in each class/stage at the start of your model. you should play around with this to see how it effects the outcome. My model is insensitive to realistic changes in this vector. Here are the “options”, the number in each class observed annually

n_options <- ddply(andre, c("YEAR"), function(df) return(table(df$CLASS)))
n_options
##    YEAR A1 A2 A3 A4 dead   J
## 1  1994  6  3  2  2    0  13
## 2  1995 31 17 13 11    5  81
## 3  1996 49 15 14  6   52  16
## 4  1997 36 19 10  4   30   2
## 5  1998 10 22 13 17    9   0
## 6  1999  8 16 24 11    3   0
## 7  2000  6 17 17  8   11   0
## 8  2001  1 11 16 14    6   0
## 9  2002 13 14  6  0    9   0
## 10 2003  4  8  6  0   14   1
## 11 2004  0  0  0  0   19   0
## 12 2011 22 23 12  6    0 173
## 13 2012  6  1  0  1  172  34

picked starting population vector from 1995, the first year with 9 observed populations

n95 <- c(81, 31, 17, 13, 11)
n = n95

you can also see that I only have J class individuals for a subset of the years. Check your data for missing bits like this. Perhaps your collaborator is forgetting to tell you that they didn’t survey for new juveniles from 1996-2003. It happens!

3. GENERATE A TRANSITION MATRIX

this matrix links each individual to its fate in the next year/cycle/season this is why you need each individual to be “dead” for a year, but no longer. If your raw data is like mine and only taken on live plants, “dead” might be something you have to add. make columns for year2, fate, and seeds2 for the whole census

trans <- subset(merge(andre, andre, by = "PLANT_UNQ", sort = FALSE), YEAR.x == 
    YEAR.y - 1)
head2(trans)
##      PLANT_UNQ YEAR.x CLASS.x SEEDS.x YEAR.y CLASS.y SEEDS.y
## 3          A_1   1994      A1   71.75   1995      A2   258.3
## 8          A_1   1995      A2  258.30   1996    dead     0.0
## 13         A_3   1996      A1   71.75   1997      A2   258.3
## .            .      .       .       .      .       .       .
## 6117     9_317   2011       J   71.75   2012       J     0.0

rename rows and columns to improve clarity (I use the names used by popbio, which are similar to Morris and Doak) I think I worked out somewhat painfully that you need these particular column names.

rownames(trans) <- 1:nrow(trans)
colnames(trans)[1:7] <- c("plant", "year", "stage", "seeds", "year2", "fate", 
    "seeds2")
head2(trans)
##     plant year stage  seeds year2 fate seeds2
## 1     A_1 1994    A1  71.75  1995   A2  258.3
## 2     A_1 1995    A2 258.30  1996 dead    0.0
## 3     A_3 1996    A1  71.75  1997   A2  258.3
## .       .    .     .      .     .    .      .
## 825 9_317 2011     J  71.75  2012    J    0.0

add individual fertility estimates from the calculations above

seedlingtrans <- 0.00305  # This is the rate at which a seed becomes a J individual (I estimated this elsewhere, see Appendix below)

adding in the number of J individuals produced by each individual

trans$J <- trans$seeds * seedlingtrans  # note that J is not an integer, which is totally fine, its a rate of J production
head2(trans)
##     plant year stage  seeds year2 fate seeds2      J
## 1     A_1 1994    A1  71.75  1995   A2  258.3 0.2188
## 2     A_1 1995    A2 258.30  1996 dead    0.0 0.7878
## 3     A_3 1996    A1  71.75  1997   A2  258.3 0.2188
## .       .    .     .      .     .    .      .      .
## 825 9_317 2011     J  71.75  2012    J    0.0 0.2188

4. GENERATE ANNUAL MATRICES — THE SIMPLE WAY for 3 easy years

NAME STAGES

stages <- c("J", "A1", "A2", "A3", "A4")

you must have a vector named stages in this way for your own clases or stages used in the “stage” vector of your data frame

SET ITERATIONS

it <- 100  # set the number of time steps for a deterministic model

Make a demographic projection matrix for each year like so: 1994

trans94 <- subset(trans, year == 1994, c(plant, stage, fate, J))
(proj94 <- projection.matrix(trans94, stage, fate, J, sort = stages))  #this gives you a projection matrix for 1994
##     
##            J     A1     A2     A3    A4
##   J  0.07692 0.1459 0.7878 0.7878 4.355
##   A1 0.53846 0.3333 0.3333 0.0000 0.000
##   A2 0.07692 0.3333 0.0000 0.0000 0.000
##   A3 0.00000 0.1667 0.3333 0.5000 0.000
##   A4 0.00000 0.0000 0.3333 0.5000 1.000

you can do a simple deterministic projection of the matrix for just this year

(p94 <- pop.projection(proj94, n, it))
## $lambda
## [1] 1.496
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.48307 0.24671 0.07983 0.06803 0.12235 
## 
## $stage.vectors
##     0     1      2      3      4      5     6      7      8      9     10
## J  81 82.29 142.62 223.46 338.78 502.17 745.4 1112.5 1665.3 2493.0 3729.6
## A1 31 59.62  69.70 108.76 167.98 256.23 383.2  570.4  850.9 1272.2 1903.7
## A2 17 16.56  26.20  34.21  53.44  82.05 124.0  185.1  275.7  411.7  615.8
## A3 13 17.33  24.12  32.41  45.74  68.68 104.4  157.4  235.5  351.4  525.0
## A4 11 23.17  37.35  58.15  85.76 126.44 188.1  281.7  422.1  631.7  944.7
##        11   12    13    14    15    16    17    18     19     20     21
## J  5577.3 8340 12472 18653 27896 41720 62394 93313 139553 208707 312129
## A1 2848.1 4260  6370  9526 14247 21307 31865 47656  71271 106588 159407
## A2  921.5 1378  2061  3083  4610  6895 10311 15421  23063  34492  51584
## A3  785.1 1174  1757  2627  3929  5876  8787 13142  19654  29393  43959
## A4 1412.4 2112  3159  4724  7065 10567 15803 23633  35345  52859  79053
##        22     23      24      25      26      27      28      29       30
## J  466801 698120 1044066 1561441 2335197 3492380 5222992 7811190 11681943
## A1 238400 356536  533214  797443 1192607 1783591 2667431 3989249  5966079
## A2  77146 115374  172547  258051  385925  577166  863175 1290912  1930610
## A3  65742  98319  147041  219905  328877  491848  735578 1100086  1645222
## A4 118227 176813  264431  395467  591436  884516 1322829 1978343  2958690
##          31       32       33       34       35        36        37
## J  17470807 26128282 39075876 58439515 87398601 130708056 195479056
## A1  8922507 13343961 19956421 29845615 44635295  66753845  99833010
## A2  2887304  4318077  6457855  9657977 14443886  21601401  32305748
## A3  2460494  3679766  5503236  8230306 12308748  18408218  27530217
## A4  4424838  6617520  9896762 14800999 22135477  33104480  49509056
##           38        39        40        41        42        43        44
## J  292346643 437215941 653873694 977893912 1.462e+09 2.187e+09 3.271e+09
## A1 149304206 223290332 333939503 499419705 7.469e+08 1.117e+09 1.671e+09
## A2  48314520  72256272 108062106 161611144 2.417e+08 3.615e+08 5.406e+08
## A3  41172526  61575137  92088048 137721310 2.060e+08 3.080e+08 4.607e+08
## A4  74042747 110733850 165606843 247671569 3.704e+08 5.540e+08 8.285e+08
##           45        46        47        48        49        50        51
## J  4.892e+09 7.316e+09 1.094e+10 1.636e+10 2.447e+10 3.660e+10 5.474e+10
## A1 2.498e+09 3.736e+09 5.588e+09 8.357e+09 1.250e+10 1.869e+10 2.795e+10
## A2 8.085e+08 1.209e+09 1.808e+09 2.704e+09 4.044e+09 6.049e+09 9.046e+09
## A3 6.890e+08 1.030e+09 1.541e+09 2.305e+09 3.447e+09 5.154e+09 7.709e+09
## A4 1.239e+09 1.853e+09 2.771e+09 4.144e+09 6.198e+09 9.270e+09 1.386e+10
##           52        53        54        55        56        57        58
## J  8.186e+10 1.224e+11 1.831e+11 2.738e+11 4.095e+11 6.124e+11 9.159e+11
## A1 4.181e+10 6.252e+10 9.351e+10 1.398e+11 2.091e+11 3.128e+11 4.678e+11
## A2 1.353e+10 2.023e+10 3.026e+10 4.525e+10 6.768e+10 1.012e+11 1.514e+11
## A3 1.153e+10 1.724e+10 2.579e+10 3.856e+10 5.767e+10 8.625e+10 1.290e+11
## A4 2.073e+10 3.101e+10 4.637e+10 6.935e+10 1.037e+11 1.551e+11 2.320e+11
##           59        60        61        62        63        64        65
## J  1.370e+12 2.049e+12 3.064e+12 4.582e+12 6.852e+12 1.025e+13 1.533e+13
## A1 6.996e+11 1.046e+12 1.565e+12 2.340e+12 3.500e+12 5.234e+12 7.827e+12
## A2 2.264e+11 3.386e+11 5.063e+11 7.572e+11 1.132e+12 1.694e+12 2.533e+12
## A3 1.929e+11 2.885e+11 4.315e+11 6.453e+11 9.651e+11 1.443e+12 2.158e+12
## A4 3.469e+11 5.188e+11 7.759e+11 1.160e+12 1.736e+12 2.596e+12 3.882e+12
##           66        67        68        69        70        71        72
## J  2.292e+13 3.428e+13 5.127e+13 7.667e+13 1.147e+14 1.715e+14 2.565e+14
## A1 1.171e+13 1.751e+13 2.618e+13 3.916e+13 5.856e+13 8.758e+13 1.310e+14
## A2 3.788e+12 5.665e+12 8.473e+12 1.267e+13 1.895e+13 2.834e+13 4.238e+13
## A3 3.228e+12 4.828e+12 7.220e+12 1.080e+13 1.615e+13 2.415e+13 3.612e+13
## A4 5.805e+12 8.682e+12 1.298e+13 1.942e+13 2.904e+13 4.343e+13 6.495e+13
##           73        74        75        76        77        78        79
## J  3.836e+14 5.736e+14 8.579e+14 1.283e+15 1.919e+15 2.870e+15 4.292e+15
## A1 1.959e+14 2.930e+14 4.381e+14 6.552e+14 9.799e+14 1.466e+15 2.192e+15
## A2 6.339e+13 9.480e+13 1.418e+14 2.120e+14 3.171e+14 4.742e+14 7.092e+14
## A3 5.402e+13 8.079e+13 1.208e+14 1.807e+14 2.702e+14 4.041e+14 6.044e+14
## A4 9.714e+13 1.453e+14 2.173e+14 3.249e+14 4.860e+14 7.268e+14 1.087e+15
##           80        81        82        83        84        85        86
## J  6.418e+15 9.599e+15 1.436e+16 2.147e+16 3.211e+16 4.802e+16 7.181e+16
## A1 3.278e+15 4.902e+15 7.331e+15 1.096e+16 1.640e+16 2.452e+16 3.668e+16
## A2 1.061e+15 1.586e+15 2.372e+15 3.548e+15 5.306e+15 7.936e+15 1.187e+16
## A3 9.039e+14 1.352e+15 2.022e+15 3.024e+15 4.522e+15 6.763e+15 1.011e+16
## A4 1.626e+15 2.431e+15 3.636e+15 5.437e+15 8.132e+15 1.216e+16 1.819e+16
##           87        88        89        90        91        92        93
## J  1.074e+17 1.606e+17 2.402e+17 3.592e+17 5.373e+17 8.035e+17 1.202e+18
## A1 5.485e+16 8.203e+16 1.227e+17 1.835e+17 2.744e+17 4.104e+17 6.137e+17
## A2 1.775e+16 2.654e+16 3.970e+16 5.937e+16 8.879e+16 1.328e+17 1.986e+17
## A3 1.513e+16 2.262e+16 3.383e+16 5.059e+16 7.566e+16 1.132e+17 1.692e+17
## A4 2.720e+16 4.068e+16 6.084e+16 9.099e+16 1.361e+17 2.035e+17 3.043e+17
##           94        95        96        97        98        99
## J  1.797e+18 2.688e+18 4.020e+18 6.011e+18 8.990e+18 1.345e+19
## A1 9.178e+17 1.373e+18 2.053e+18 3.070e+18 4.591e+18 6.867e+18
## A2 2.970e+17 4.442e+17 6.643e+17 9.935e+17 1.486e+18 2.222e+18
## A3 2.531e+17 3.785e+17 5.661e+17 8.466e+17 1.266e+18 1.894e+18
## A4 4.552e+17 6.807e+17 1.018e+18 1.522e+18 2.277e+18 3.405e+18
## 
## $pop.sizes
##   [1] 1.530e+02 1.990e+02 3.000e+02 4.570e+02 6.917e+02 1.036e+03 1.545e+03
##   [8] 2.307e+03 3.449e+03 5.160e+03 7.719e+03 1.154e+04 1.726e+04 2.582e+04
##  [15] 3.861e+04 5.775e+04 8.636e+04 1.292e+05 1.932e+05 2.889e+05 4.320e+05
##  [22] 6.461e+05 9.663e+05 1.445e+06 2.161e+06 3.232e+06 4.834e+06 7.230e+06
##  [29] 1.081e+07 1.617e+07 2.418e+07 3.617e+07 5.409e+07 8.089e+07 1.210e+08
##  [36] 1.809e+08 2.706e+08 4.047e+08 6.052e+08 9.051e+08 1.354e+09 2.024e+09
##  [43] 3.027e+09 4.528e+09 6.771e+09 1.013e+10 1.514e+10 2.265e+10 3.387e+10
##  [50] 5.066e+10 7.576e+10 1.133e+11 1.695e+11 2.534e+11 3.790e+11 5.668e+11
##  [57] 8.477e+11 1.268e+12 1.896e+12 2.836e+12 4.241e+12 6.342e+12 9.485e+12
##  [64] 1.419e+13 2.121e+13 3.173e+13 4.745e+13 7.096e+13 1.061e+14 1.587e+14
##  [71] 2.374e+14 3.550e+14 5.309e+14 7.940e+14 1.187e+15 1.776e+15 2.656e+15
##  [78] 3.972e+15 5.940e+15 8.884e+15 1.329e+16 1.987e+16 2.972e+16 4.444e+16
##  [85] 6.646e+16 9.940e+16 1.487e+17 2.223e+17 3.325e+17 4.973e+17 7.437e+17
##  [92] 1.112e+18 1.663e+18 2.488e+18 3.720e+18 5.564e+18 8.321e+18 1.244e+19
##  [99] 1.861e+19 2.783e+19
## 
## $pop.changes
##  [1] 1.300 1.508 1.523 1.514 1.497 1.492 1.493 1.495 1.496 1.496 1.496
## [12] 1.495 1.495 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [23] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [34] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [45] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [56] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [67] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [78] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [89] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
(l94 <- p94$lambda)  # wow! if we looked only at 1994 based on these estimates the population would be booming!
## [1] 1.496

stable.stage shows the proportion of the population in each stage class at the mythical equilibrium, 48% of plants are juveniles in 100 years note that the population has gone from around 200 individuals in 1994 to 2.8e+19 in 100 years, nice!

Now make some matrices for other years 1995

trans95 <- subset(trans, year == 1995, c(plant, stage, fate, J))
(proj95 <- projection.matrix(trans95, stage, fate, J, sort = stages))
##     
##           J     A1     A2      A3      A4
##   J  0.2000 0.2626 1.0298 3.52160 6.21499
##   A1 0.2125 0.7419 0.4706 0.07692 0.00000
##   A2 0.0000 0.1290 0.4706 0.15385 0.09091
##   A3 0.0000 0.0000 0.0000 0.76923 0.36364
##   A4 0.0000 0.0000 0.0000 0.00000 0.54545

2011

trans11 <- subset(trans, year == 2011, c(plant, stage, fate, J))
(proj11 <- projection.matrix(trans11, stage, fate, J, sort = stages))
##     
##           J      A1     A2     A3    A4
##   J  0.1573 0.98329 4.5710 8.0648 7.922
##   A1 0.0000 0.05556 0.3077 0.1429 0.000
##   A2 0.0000 0.00000 0.0000 0.0000 0.500
##   A3 0.0000 0.00000 0.0000 0.0000 0.000
##   A4 0.0000 0.00000 0.0000 0.1429 0.000

p95 <- pop.projection(proj95, n, it)
(l95 <- p95$lambda)  # lambda is much lower in 1995
## [1] 0.9954

(p11 <- pop.projection(proj11, n, pi))
## $lambda
## [1] 0.3065
## 
## $stable.stage
##        J       A1       A2       A3       A4 
## 0.969161 0.021632 0.009207 0.000000 0.000000 
## 
## $stage.vectors
##     0       1       2
## J  81 312.916 97.7480
## A1 31   8.810  2.1818
## A2 17   5.500  0.9286
## A3 13   0.000  0.0000
## A4 11   1.857  0.0000
## 
## $pop.sizes
## [1] 153.0 329.1 100.9
## 
## $pop.changes
## [1] 2.1509 0.3065
(l11 <- p11$lambda)  # and based on 2011 alone extinction is eminent. The gist here is we need lots of years of data to make any decent estimation
## [1] 0.3065
# of what the population is really likely to do (ie more than the three
# here)

5. DETERMINISTIC ANALAYSIS OF THESE 3 WELL BEHAVED YEARS - basic non-stochastic PVA

thesearethemeanprojmats <- list(proj94, proj95, proj11)  # make a list of the three matrices
(meanxprojmat <- mean(thesearethemeanprojmats))  # make a mean of the three projection matrices for deterministic analysis
##     
##            J      A1     A2      A3     A4
##   J  0.14475 0.46393 2.1295 4.12473 6.1639
##   A1 0.25032 0.37694 0.3705 0.07326 0.0000
##   A2 0.02564 0.15412 0.1569 0.05128 0.1970
##   A3 0.00000 0.05556 0.1111 0.42308 0.1212
##   A4 0.00000 0.00000 0.1111 0.21429 0.5152

n  # recall what n is, our starting population vector, ie the # of individuals in each class at the start of the projection
## [1] 81 31 17 13 11
(pprojme <- pop.projection(meanxprojmat, n))  # do the deterministic projection, lambda is the dominant left eigenvalue
## $lambda
## [1] 1.089
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.61436 0.25469 0.06677 0.03730 0.02688 
## 
## $stage.vectors
##     0      1       2       3       4       5      6      7      8      9
## J  81 183.73 177.921 182.074 196.519 214.113 233.26 253.98 276.50 301.01
## A1 31  39.21  66.116  75.791  82.085  88.980  96.75 105.30 114.63 124.79
## A2 17  12.35  15.265  19.380  21.551  23.374  25.38  27.61  30.05  32.71
## A3 13  10.44   9.224  10.355  11.748  12.972  14.16  15.42  16.79  18.28
## A4 11  10.34   8.938   8.277   8.636   9.361  10.20  11.11  12.10  13.17
##        10     11     12     13     14     15     16     17     18     19
## J  327.68 356.72 388.33 422.74 460.21 500.99 545.38 593.72 646.33 703.61
## A1 135.85 147.88 160.99 175.26 190.79 207.69 226.10 246.14 267.95 291.69
## A2  35.61  38.77  42.20  45.94  50.02  54.45  59.27  64.53  70.24  76.47
## A3  19.90  21.66  23.58  25.67  27.94  30.42  33.11  36.05  39.24  42.72
## A4  14.33  15.60  16.99  18.49  20.13  21.92  23.86  25.97  28.27  30.78
## 
## $pop.sizes
##  [1]  153.0  256.1  277.5  295.9  320.5  348.8  379.7  413.4  450.1  490.0
## [11]  533.4  580.6  632.1  688.1  749.1  815.5  887.7  966.4 1052.0 1145.3
## 
## $pop.changes
##  [1] 1.674 1.083 1.066 1.083 1.088 1.089 1.089 1.089 1.089 1.089 1.089
## [12] 1.089 1.089 1.089 1.089 1.089 1.089 1.089 1.089
(DetLamb <- pprojme$lambda)
## [1] 1.089

6. RUN STOCHASTIC ANALYSES ON THESE 3 WELL BEHAVED YEARS

For a stochastic analysis, include all of the annual matrices, then make a random draw with replacement

For a series of time step, or until stable stage distribution is reached.

thesearethemeanprojmats  # our list of projection matrices
## [[1]]
##     
##            J     A1     A2     A3    A4
##   J  0.07692 0.1459 0.7878 0.7878 4.355
##   A1 0.53846 0.3333 0.3333 0.0000 0.000
##   A2 0.07692 0.3333 0.0000 0.0000 0.000
##   A3 0.00000 0.1667 0.3333 0.5000 0.000
##   A4 0.00000 0.0000 0.3333 0.5000 1.000
## 
## [[2]]
##     
##           J     A1     A2      A3      A4
##   J  0.2000 0.2626 1.0298 3.52160 6.21499
##   A1 0.2125 0.7419 0.4706 0.07692 0.00000
##   A2 0.0000 0.1290 0.4706 0.15385 0.09091
##   A3 0.0000 0.0000 0.0000 0.76923 0.36364
##   A4 0.0000 0.0000 0.0000 0.00000 0.54545
## 
## [[3]]
##     
##           J      A1     A2     A3    A4
##   J  0.1573 0.98329 4.5710 8.0648 7.922
##   A1 0.0000 0.05556 0.3077 0.1429 0.000
##   A2 0.0000 0.00000 0.0000 0.0000 0.500
##   A3 0.0000 0.00000 0.0000 0.0000 0.000
##   A4 0.0000 0.00000 0.0000 0.1429 0.000

stochme <- stoch.growth.rate(thesearethemeanprojmats, prob = NULL, maxt = 50000, 
    verbose = TRUE)
## [1] Calculating stochastic growth at time 1
## [1] Calculating stochastic growth at time 10000
## [1] Calculating stochastic growth at time 20000
## [1] Calculating stochastic growth at time 30000
## [1] Calculating stochastic growth at time 40000
## [1] Calculating stochastic growth at time 50000

note that these stochastic approximations of lambda are in log form (not immediately comparable to pop.project$lambda)

exp(stochme$approx)  # is the analytic approximation of lambda via Tuljapakar's method
## [1] 0.946
stochme$approx  # this is more accurate (perhaps) when there is a lot of covariation in matrix elements
## [1] -0.05556
exp(stochme$sim)  # gives stochastic growth rate by simulation, random draws of whole
## [1] 0.8884
stochme$sim
## [1] -0.1183

Tuljapurkar’s approximation takes into account how stochastic variation in the matrix elements affects the long-term stochastic growth rate (Caswell 2001). It can be more accurate in cases where there is covariation between matrix elements within the same year but may not be as accurate when there is a high level of temporal variation (Morris and Doak 2002, Stubben et al. 2012).

Fun with stochastic analyses!

Its easy to give years “weights” in the stochastic model. For example, you can increase the drought rate by weighting drought years (2011)

yearweight <- c(1, 1, 2)
moredrought <- stoch.projection(thesearethemeanprojmats, n, tmax = 50, prob = yearweight, 
    nreps = 500)

the output is population sizes, which are fun to graph when comparing models

yearweight <- c(1, 1, 0)
nodrought <- stoch.projection(thesearethemeanprojmats, n, tmax = 50, prob = yearweight, 
    nreps = 500)

par(mfrow = c(2, 1))
hist(log(apply(moredrought, 1, sum)), col = "blue", density = 50, ylim = c(0, 
    150), xlim = c(-1.3, 25), xlab = "", main = "More drought")
abline(v = log10(200), lty = 3)  # puts a line at the starting population size for reference
hist(log(apply(nodrought, 1, sum)), col = "green3", density = 50, ylim = c(0, 
    150), xlim = c(-1.3, 25), xlab = "", main = "No drought")
abline(v = log10(200), lty = 3)

you can get fancy and put these on the same graph too to compare outcomes. y axis is frequency of final population size at tmax.

7. QUASI EXTINCTION BASED ON THESE 3 WELL BEHAVED YEARS

Another useful way to think about populations. Since our ability to really estimate lambda is based on the assumption of equilibrium at stable stage, it might be more realistic to think about comparing extinction probabilities for different scenarios that comparing lambdas. these are based on stochastic runs

obsd <- stoch.quasi.ext(thesearethemeanprojmats, n, prob = c(1, 1, 1), Nx = 10, 
    tmax = 50, maxruns = 10, nreps = 500, sumweight = c(1, 1, 1), verbose = FALSE)

drt <- stoch.quasi.ext(thesearethemeanprojmats, n, prob = c(1, 1, 2), Nx = 10, 
    tmax = 50, maxruns = 10, nreps = 500, sumweight = c(1, 1, 1), verbose = FALSE)

par(mfrow = c(2, 1))
matplot(obsd, ylab = "Quasi-extinction probability", ylim = c(0, 1.1), type = "l", 
    lty = 1, col = rainbow(10), las = 1, main = "Observed climate", xlab = "Years")
matplot(drt, ylab = "Quasi-extinction probability", ylim = c(0, 1.1), type = "l", 
    lty = 1, col = rainbow(10), las = 1, main = "Double drought", xlab = "Years")

8. SENSITIVITY AND ELASTICITY

SENSITIVITY is a measure of the amount of change is lambda give a small change in a matrix element.

ELASTICITY is a measure of ``proportional’’ effect, i.e., the effect that a change in a given matrix element has as a proportional to the change in that element

meanxprojmat
##     
##            J      A1     A2      A3     A4
##   J  0.14475 0.46393 2.1295 4.12473 6.1639
##   A1 0.25032 0.37694 0.3705 0.07326 0.0000
##   A2 0.02564 0.15412 0.1569 0.05128 0.1970
##   A3 0.00000 0.05556 0.1111 0.42308 0.1212
##   A4 0.00000 0.00000 0.1111 0.21429 0.5152

for an overall look at sensitivity and elasticity use the mean projection matrix you could do separate analyses by year or type of year too to examine how sensitivity and elasticity vary among years

(eigout <- eigen.analysis(meanxprojmat))  # do the associated sensitivity analysis
## $lambda1
## [1] 1.089
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.61436 0.25469 0.06677 0.03730 0.02688 
## 
## $sensitivities
##     
##           J      A1      A2      A3       A4
##   J  0.2255 0.09351 0.02451 0.01369 0.009867
##   A1 0.6930 0.28729 0.07532 0.04208 0.000000
##   A2 1.5373 0.63734 0.16708 0.09334 0.067253
##   A3 0.0000 1.13131 0.29658 0.16569 0.119377
##   A4 0.0000 0.00000 0.38356 0.21428 0.154386
## 
## $elasticities
##     
##            J      A1      A2       A3      A4
##   J  0.02999 0.03985 0.04795 0.051889 0.05587
##   A1 0.15935 0.09948 0.02564 0.002832 0.00000
##   A2 0.03621 0.09023 0.02408 0.004397 0.01217
##   A3 0.00000 0.05773 0.03027 0.064394 0.01329
##   A4 0.00000 0.00000 0.03915 0.042180 0.07306
## 
## $repro.value
##      J     A1     A2     A3     A4 
##  1.000  3.072  6.816 12.099 15.647 
## 
## $damping.ratio
## [1] 4.215

colSums(eigout$sensitivities)  # this gives the cumulative sensitivity of each stage/class, fun to graph
##      J     A1     A2     A3     A4 
## 2.4559 2.1494 0.9471 0.5291 0.3509

calculate fertility and survival sums

(projsums <- colSums(meanxprojmat))
##      J     A1     A2     A3     A4 
## 0.4207 1.0505 2.8792 4.8866 6.9973
(fert_row <- meanxprojmat[1, ])
##      J     A1     A2     A3     A4 
## 0.1448 0.4639 2.1295 4.1247 6.1639
(surv_row <- projsums - fert_row)
##      J     A1     A2     A3     A4 
## 0.2760 0.5866 0.7496 0.7619 0.8333

can do this for sens and elas too to make a bar chart

make x where columes are name, sens and elas

for my data I’m summing elasticity and sensitivity for each class. You could also pick the vital rates that are most meaningful in your own analysis

(eigout <- eigen.analysis(meanxprojmat))
## $lambda1
## [1] 1.089
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.61436 0.25469 0.06677 0.03730 0.02688 
## 
## $sensitivities
##     
##           J      A1      A2      A3       A4
##   J  0.2255 0.09351 0.02451 0.01369 0.009867
##   A1 0.6930 0.28729 0.07532 0.04208 0.000000
##   A2 1.5373 0.63734 0.16708 0.09334 0.067253
##   A3 0.0000 1.13131 0.29658 0.16569 0.119377
##   A4 0.0000 0.00000 0.38356 0.21428 0.154386
## 
## $elasticities
##     
##            J      A1      A2       A3      A4
##   J  0.02999 0.03985 0.04795 0.051889 0.05587
##   A1 0.15935 0.09948 0.02564 0.002832 0.00000
##   A2 0.03621 0.09023 0.02408 0.004397 0.01217
##   A3 0.00000 0.05773 0.03027 0.064394 0.01329
##   A4 0.00000 0.00000 0.03915 0.042180 0.07306
## 
## $repro.value
##      J     A1     A2     A3     A4 
##  1.000  3.072  6.816 12.099 15.647 
## 
## $damping.ratio
## [1] 4.215
eigout$sensitivities
##     
##           J      A1      A2      A3       A4
##   J  0.2255 0.09351 0.02451 0.01369 0.009867
##   A1 0.6930 0.28729 0.07532 0.04208 0.000000
##   A2 1.5373 0.63734 0.16708 0.09334 0.067253
##   A3 0.0000 1.13131 0.29658 0.16569 0.119377
##   A4 0.0000 0.00000 0.38356 0.21428 0.154386

(fert_row_s <- eigout$sensitivities[1, ])
##        J       A1       A2       A3       A4 
## 0.225549 0.093506 0.024513 0.013695 0.009867
(surv_row_s <- eigout$sensitivities[2, ] + eigout$sensitivities[3, ] + eigout$sensitivities[4, 
    ] + eigout$sensitivities[5, ])
##      J     A1     A2     A3     A4 
## 2.2303 2.0559 0.9225 0.5154 0.3410
sensme <- t(cbind(fert_row_s, surv_row_s))


(fert_row_e <- eigout$elasticities[1, ])
##       J      A1      A2      A3      A4 
## 0.02999 0.03985 0.04795 0.05189 0.05587
(surv_row_e <- eigout$elasticities[2, ] + eigout$elasticities[3, ] + eigout$elasticities[4, 
    ] + eigout$elasticities[5, ])
##       J      A1      A2      A3      A4 
## 0.19556 0.24744 0.11913 0.11380 0.09852
(survtable <- t(rbind(surv_row_s, surv_row_e)))
##    surv_row_s surv_row_e
## J      2.2303    0.19556
## A1     2.0559    0.24744
## A2     0.9225    0.11913
## A3     0.5154    0.11380
## A4     0.3410    0.09852
(ferttable <- t(rbind(fert_row_s, fert_row_e)))
##    fert_row_s fert_row_e
## J    0.225549    0.02999
## A1   0.093506    0.03985
## A2   0.024513    0.04795
## A3   0.013695    0.05189
## A4   0.009867    0.05587

par(mfrow = c(2, 2))
mynames <- c("Sensitivity", "Elasticity")
barplot(t(survtable[, 1:2]), beside = TRUE, las = 1, ylim = c(0, 3.5), xlab = "Stage class", 
    main = "Growth and survival")
abline(h = 0)
barplot(t(ferttable[, 1:2]), beside = TRUE, las = 1, ylim = c(0, 0.25), xlab = "Stage class", 
    main = "Fertility")
abline(h = 0)
legend("topright", mynames, fill = grey.colors(2))

APPENDIX A. HOW I GENERATED A FECUNDITY ESTIMATE

I’m including this for plant folks who might like to see how I made fecundity estimates from real data I also have developed a script to simulate juvenile numbers and transition rates for the years in my dataset that are missing these data, and bootstraps of the whole model.

Load up the raw-ish data

andre <- read.csv("D__composite9_13_2012.csv")
str(andre)
## 'data.frame':    1693 obs. of  15 variables:
##  $ PLANT_UNQ   : Factor w/ 619 levels "1_701","1_702",..: 462 464 473 477 482 487 474 479 483 463 ...
##  $ YEAR        : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS       : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ OBS_ID      : int  1 3 20 24 11 16 21 26 12 2 ...
##  $ COHORT      : Factor w/ 18 levels "1","2","3","4",..: 10 10 11 11 11 11 11 11 11 10 ...
##  $ PLANT.NUM   : Factor w/ 460 levels "1","10","11",..: 1 23 5 9 41 168 6 11 42 12 ...
##  $ xDIAM_cm    : num  12 11 17 12 18 12 23 22 29 32 ...
##  $ SEEDS_DC    : num  26.7 26.7 26.7 26.7 26.7 ...
##  $ MEDIAN_INFL : int  5 0 5 0 5 5 18 18 18 18 ...
##  $ INFL_CLASS  : int  1 0 1 0 1 1 2 2 2 2 ...
##  $ MAX_INFL    : int  111 0 111 0 111 111 123 123 123 123 ...
##  $ OBS_INLF    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ SEEDS_YR_INF: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ CAGED       : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
##  $ WATERED     : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...

BEST ESTIMATE OF SEED PRODUCTION: average of seeds/fruit, weighed average of fruits/infl where 0.14 is 1/8 drought years

andre$SEEDS <- andre$MEDIAN_INFL * 14.35  # = 14.35 seeds/infl

The big crux Part 2: How many seeds makes a J plant the next year? ie what is the transition rate or fecundity rate? observed juveniles in each year (sadly not all years have the same # of cohorts, so I adjust for that below)

seedling_yr <- ddply(andre, c("YEAR"), function(df) return(c(sdlgs = sum(df$CLASS == 
    "J"))))  # shows the OBSERVED # of J plants each year,

and some of these are not observation years and are omitted below

seedling_yr
##    YEAR sdlgs
## 1  1994    13
## 2  1995    81
## 3  1996    16
## 4  1997     2
## 5  1998     0
## 6  1999     0
## 7  2000     0
## 8  2001     0
## 9  2002     0
## 10 2003     1
## 11 2004     0
## 12 2011   362
## 13 2012    57

seedlings for 9 cohorts based on census from each year

sdl94 = seedling_yr$sdlgs[(seedling_yr$YEAR == "1994")] * 9/2  # this year had only 2 cohorts
sdl95 = seedling_yr$sdlgs[(seedling_yr$YEAR == "1995")]
sdl11 = seedling_yr$sdlgs[(seedling_yr$YEAR == "2011")]
sdl12 = 0  # in 2012 observed seedlings were 0

seedling_pick = c(sdl94, sdl95, sdl11, sdl12)
seedling_pick
## [1]  58.5  81.0 362.0   0.0

estimate seedlings as the mean of the 4 OBSERVED years: in the other years seedlings where not surveyed for

(seedlings = mean(seedling_pick))  # = 125.375 seedlings/year on average
## [1] 125.4

what’s the annual total seed production rate? add up all the seeds estimated to be produced in each year

seed_yr <- ddply(andre, c("YEAR"), function(df) return(c(sumseeds = sum(df$SEEDS))))
seed_yr
##    YEAR sumseeds
## 1  1994     4434
## 2  1995    45834
## 3  1996    13302
## 4  1997    15426
## 5  1998    51976
## 6  1999    45030
## 7  2000    18913
## 8  2001    45576
## 9  2002     2597
## 10 2003     6257
## 11 2004        0
## 12 2011   158065
## 13 2012     3071

adjust so that 1994 has an estimate of all cohorts based on the observed 2 cohorts

seed_yr$sumseeds[(seedling_yr$YEAR == "1994")] <- (seed_yr$sumseeds[(seedling_yr$YEAR == 
    "1994")] * 9/2)

get mean seeds/year

avg_seeds_p_yr <- mean(seed_yr$sumseeds)

for each of the 4 years in which seedlings were observed, calculate an estimate of the transition rate from seed –> J

str(sdl.trans94 <- sdl94/avg_seeds_p_yr)
##  num 0.00179
str(sdl.trans95 <- sdl95/avg_seeds_p_yr)
##  num 0.00247
str(sdl.trans11 <- sdl11/avg_seeds_p_yr)
##  num 0.011
sdl.trans12 <- 0  # no seedlings observed this year

looks at the options for transition rate

str(seedlingtrans_pick <- c(sdl.trans94, sdl.trans95, sdl.trans11, sdl.trans12))
##  num [1:4] 0.00179 0.00247 0.01105 0

pick the mean for this analysis

(seedlingtrans <- mean(seedlingtrans_pick))  # mean seedling transition rate  0.00305
## [1] 0.003826

Now remove watered and caged plants from main dataset, these have different survival and transition rates. I included them thus far because we needed to get seed production estimates for them. Since they might have contributed to the observed juveniles

andre <- andre[(andre$CAGED == "N"), ]
andre <- andre[(andre$WATERED == "N"), ]
head2(andre)
##      PLANT_UNQ YEAR CLASS OBS_ID COHORT PLANT.NUM xDIAM_cm SEEDS_DC
## 1          A_1 1994    A1      1      A         1     12.0     26.7
## 2          A_3 1994    A1      3      A         3     11.0     26.7
## 3         B_13 1994    A1     20      B        13     17.0     26.7
## .            .    .     .      .      .         .        .        .
## 1675    6_772A 2012     J    608      6      772A      4.7      0.0
##      MEDIAN_INFL INFL_CLASS MAX_INFL OBS_INLF SEEDS_YR_INF CAGED WATERED
## 1              5          1      111       NA           NA     N       N
## 2              0          0        0       NA           NA     N       N
## 3              5          1      111       NA           NA     N       N
## .              .          .        .        .            .     .       .
## 1675           0          0        0        0            0     N       N
##      SEEDS
## 1    71.75
## 2     0.00
## 3    71.75
## .        .
## 1675  0.00

reduce datafile to include only PLANT_UNQ, YEAR, CLASS and SEEDS

str(andre)
## 'data.frame':    1221 obs. of  16 variables:
##  $ PLANT_UNQ   : Factor w/ 619 levels "1_701","1_702",..: 462 464 473 477 482 487 474 479 483 463 ...
##  $ YEAR        : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS       : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ OBS_ID      : int  1 3 20 24 11 16 21 26 12 2 ...
##  $ COHORT      : Factor w/ 18 levels "1","2","3","4",..: 10 10 11 11 11 11 11 11 11 10 ...
##  $ PLANT.NUM   : Factor w/ 460 levels "1","10","11",..: 1 23 5 9 41 168 6 11 42 12 ...
##  $ xDIAM_cm    : num  12 11 17 12 18 12 23 22 29 32 ...
##  $ SEEDS_DC    : num  26.7 26.7 26.7 26.7 26.7 ...
##  $ MEDIAN_INFL : int  5 0 5 0 5 5 18 18 18 18 ...
##  $ INFL_CLASS  : int  1 0 1 0 1 1 2 2 2 2 ...
##  $ MAX_INFL    : int  111 0 111 0 111 111 123 123 123 123 ...
##  $ OBS_INLF    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ SEEDS_YR_INF: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ CAGED       : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
##  $ WATERED     : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
##  $ SEEDS       : num  71.8 0 71.8 0 71.8 ...
str(andre[, c(1:3, 16)])
## 'data.frame':    1221 obs. of  4 variables:
##  $ PLANT_UNQ: Factor w/ 619 levels "1_701","1_702",..: 462 464 473 477 482 487 474 479 483 463 ...
##  $ YEAR     : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS    : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ SEEDS    : num  71.8 0 71.8 0 71.8 ...

Go back to step 2.