Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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:
- Estimating the population growth rate – deterministic and stochastic methods
- Sensitivity and elasticity analysis – which transitions/rate are most variable, sensitive to change
- Projecting the population stochastically for different scenarios, observing variation in growth rate
- Quasi extinction probability
- 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.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.