vignettes/a02_VitalRates.Rmd
a02_VitalRates.RmdThe transition rates that make up matrix population models (MPMs) generally reflect products of two or more vital rates (sometimes called ‘lower-level vital rates’). Assuming a post-breeding census design, we can retroactively break apart each transition rate into at least two vital rate components: survival, and ‘something’ conditional on survival. That ‘something’ might be growth, shrinkage, stasis, dormancy, fecundity, or clonality.
The vr_ group of functions can be used to derive vital
rates from any MPM that was initially estimated using a post-breeding
census design. The vital rates may be summarized across the entire MPM
(vr_ functions), summarized within stage classes
(vr_vec_ functions), or not summarized
(vr_mat_ functions).
We’ll start by creating an MPM based on a set of underlying vital
rates. Our MPM will be separated into a growth/survival component
(matU) and a sexual reproduction component
(matF), and consist of 4 stage classes: ‘seed’, ‘small’,
‘large’, and ‘dormant’. First we’ll specify the vital rates.
## survival probabilities
s1 <- 0.20 # seed
s2 <- 0.40 # small
s3 <- 0.50 # large
s4 <- 0.30 # dormant
## matU vital rates conditional on survival (growth, shrinkage, stasis, etc.)
g12 <- 0.15 # seed to small (growth)
g13 <- 0.05 # seed to large (growth)
g11 <- 1 - g12 - g13 # seed to seed (stasis)
g23 <- 0.45 # small to large (growth)
g24 <- 0.10 # small to dormant (enter dormancy)
g22 <- 1 - g23 - g24 # small to small (stasis)
g34 <- 0.22 # large to dormant (enter dormancy)
g33 <- 1 - g34 # large to large (stasis)
g43 <- 0.50 # dormant to large (exit dormancy)
g44 <- 1 - g43 # dormant to dormant (stasis)
## matF vital rates conditional on survival (i.e. fecundity)
f2 <- 0.4 # small
f3 <- 1.1 # largeNext we’ll use the vital rates to construct an MPM, based on a post-breeding census design.
Given a post-breeding census design, stage-specific vital rates of
survival can be obtained by taking column sums of the U
matrix (matU in our code). This works out algebraically
because the other vital rates within matU, which are
conditional on survival, must sum to 1 within a given stage class. We
can obtain a vector of stage-specific vital rates of survival using the
vr_vec_survival() function, which, for simple cases, is
equivalent to colSums().
(surv <- vr_vec_survival(matU))
#> [1] 0.2 0.4 0.5 0.3
# equivalent to...
(surv <- colSums(matU))
#> [1] 0.2 0.4 0.5 0.3As expected, these values match the survival probabilities that we initially specified.
To obtain the other set of vital rates — vital rates conditional on
survival (sometimes called ‘survival-independent vital rates’) — we
simply need to divide each column of matU or
matF by the corresponding survival probability. This
procedure is implemented by the vr_mat_ functions.
# matU vital rates conditional on survival
vr_mat_U(matU)
#> [,1] [,2] [,3] [,4]
#> [1,] 0.80 NA NA NA
#> [2,] 0.15 0.45 NA NA
#> [3,] 0.05 0.45 0.78 0.5
#> [4,] NA 0.10 0.22 0.5
# matF vital rates conditional on survival
vr_mat_R(matU, matF)
#> [,1] [,2] [,3] [,4]
#> [1,] NA 0.4 1.1 NA
#> [2,] NA NA NA NA
#> [3,] NA NA NA NA
#> [4,] NA NA NA NAWe’ve now recovered all the vital rates that we initially specified when we constructed our MPM.
The vital rates produced by vr_mat_U() represent a
variety of processes, including growth, shrinkage, stasis, and dormancy.
To summarize these processes within stage classes we can use the
vr_vec_ group of functions, which simply sum the vital
rates corresponding to a given process within stage classes.
We’ll start by examining vital rates of growth, which will be found
in the bottom-left triangle of the U matrix. Given the
results from vr_mat_U() above, our baseline expectation for
the stage-specific vital rates of growth would be 0.20 (0.15 + 0.05),
0.55 (0.45 + 0.10), 0.22, and NA.
vr_vec_growth(matU)
#> [1] 0.20 0.55 0.22 NAHowever, the calculations above include transitions to the dormant
stage class as a component of ‘growth’, because the dormant stage class
is represented by the right-most column (and so transitions to the
dormant stage are in the bottom-left triangle of matU). To
exclude transitions to the dormant stage from our calculation
of growth, we can set the exclude_row argument to
4, so that any survival-independent transition to
stage 4 is excluded from the calculation of growth.
vr_vec_growth(matU, exclude_row = 4)
#> [1] 0.20 0.45 NA NANext let’s examine vital rates of shrinkage, which will occur in the
top-right triangle of matU. Based again on the results from
vr_mat_U(), our naive expectation for vital rates of
shrinkage would be NA, NA, NA,
0.5.
vr_vec_shrinkage(matU)
#> [1] NA NA NA 0.5However, in this case, the transition from the dormant stage to the
large stage is being mistakenly considered as shrinkage, simply because
it’s in the top-right triangle of matU. To exclude
transitions from the dormant stage from our calculations, we
can set the exclude_col argument to 4. Note
that we use exclude_row to exclude transitions to
a given stage, and exclude_col to exclude transitions
from a given stage.
vr_vec_shrinkage(matU, exclude_col = 4)
#> [1] NA NA NA NANext let’s calculate vital rates of stasis, which occur along the
diagonal of matU. Depending on our purpose, we may or may
not consider the dormant to dormant transition as ‘stasis’. If not, we
could exclude it using one of the exclude_ arguments, as
above. But we’ll leave it in for now.
vr_vec_stasis(matU)
#> [1] 0.80 0.45 0.78 0.50Finally, let’s calculate vital rates for entering and exiting
dormancy. For these, we’ll need to explicitly specify which stage(s) are
dormant using the dorm_stages argument.
vr_vec_dorm_enter(matU, dorm_stages = 4)
#> [1] NA 0.10 0.22 NA
vr_vec_dorm_exit(matU, dorm_stages = 4)
#> [1] NA NA NA 0.5The survival-independent vital rates of fecundity produced by
vr_mat_R() are more straightforward to summarize, as they
only reflect a single process (i.e. fecundity).
vr_vec_reproduction(matU, matF)
#> [1] NA 0.4 1.1 NAIn the example MPM we’re using here, individuals produced by
fecundity always begin life in the ‘seed’ stage. However, for some life
cycles there are multiple stages in which an offspring could begin life.
In these cases, it may be desirable to weight the different types of
offspring when summing vital rates of fecundity within a stage class
(e.g. by their reproductive value). This type of weighting can be
accomplished with the weights_row argument.
In some cases, it may be desirable to summarize vital rates across
stage classes to obtain a single mean vital rate for the entire MPM.
This can be done using the vr_ group of functions.
By default, the vr_ functions take a simple average of
the stage-specific vital rates returned by the corresponding
vr_vec_ function. Here’s an example with growth.
vr_growth(matU, exclude_row = 4)
#> [1] 0.325
# equivalent to
(vec_growth <- vr_vec_growth(matU, exclude_row = 4))
#> [1] 0.20 0.45 NA NA
mean(vec_growth, na.rm = TRUE)
#> [1] 0.325Note that the two stage classes from which there are no growth
transitions do not contribute to the MPM-average vital rate of growth.
In practice, this is because the vr_vec_ functions return
NAs rather than 0s for stages that do not exhibit the relevant
transition type. This is an important issue that we’ll return to in a
later section.
Rather than taking a simple average of the given vital rate across
stage classes, we may wish to take a weighted average across
stage classes. For instance, we may wish to weight stage classes based
on the stable distribution at equilibrium. Here’s an example of how to
do that using the weights_col argument.
# calculate the stable distribution using popdemo::eigs
library(popdemo)
#> Welcome to popdemo! This is version 1.3-0
#> Use ?popdemo for an intro, or browseVignettes('popdemo') for vignettes
#> Citation for popdemo is here: doi.org/10.1111/j.2041-210X.2012.00222.x
#> Development and legacy versions are here: github.com/iainmstott/popdemo
matA <- matU + matF
w <- popdemo::eigs(matA, what = "ss")
# calculate mean vital rate of growth weighted by the stable distribution
vr_growth(matU, exclude_row = 4, weights_col = w)
#> [1] 0.2220831
# equivalent to
pos <- !is.na(vec_growth) # possible transitions
sum(vec_growth[pos] * w[pos]) / sum(w[pos]) # weighted average
#> [1] 0.2220831By default, all vr_ functions assume that a transition
rate of 0 indicates an impossible transition within the given life cycle
(e.g. tadpoles never revert to eggs), in which case a value of
NA will be used in any relevant calculations. However, a
transition rate of 0 could alternatively indicate a transition that is
generally possible, but was simply estimated to be 0 in the relevant
population and time period.
This distinction between possible and impossible transitions can be important when calculating vital rates — particularly if we want to summarize vital rates across stages (i.e. summarize at the level of the MPM).
Let’s consider the same matU we defined above, but
imagine that it reflects a life cycle with a different set of stages:
‘seed’, ‘small’, ‘medium’, and ‘large’. Now our MPM has a single
shrinkage transition (large to medium). If we summarize the vital rate
of shrinkage across stages using vr_shrinkage(), we’ll
simply recover the single shrinkage vital rate from the large to medium
transition (i.e. because the other stages, with no shrinkage, do not
contribute to the average).
vr_vec_shrinkage(matU)
#> [1] NA NA NA 0.5
vr_shrinkage(matU)
#> [1] 0.5However, imagine now that the medium to small transition, which is 0
within matU, is in fact a possible shrinkage transition in
our life cycle of interest — its rate was simply estimated to be 0 in
the relevant sample and time period. In this case, we’ll need to
manually specify the matrix of possible transitions using the
posU argument.
(pos <- matU > 0) # possible transitions based on matU
#> [,1] [,2] [,3] [,4]
#> [1,] TRUE FALSE FALSE FALSE
#> [2,] TRUE TRUE FALSE FALSE
#> [3,] TRUE TRUE TRUE TRUE
#> [4,] FALSE TRUE TRUE TRUE
pos[2, 3] <- TRUE # set medium-to-small shrinkage tr to possible
vr_vec_shrinkage(matU, posU = pos)
#> [1] NA NA 0.0 0.5
vr_shrinkage(matU, posU = pos)
#> [1] 0.25Now, instead of returning NA for the vital rate of
shrinkage from the medium stage class, vr_vec_shrinkage()
returns a value of 0. Correspondingly, vr_shrinkage()
incorporates this 0 into its calculation and returns an MPM-average
shrinkage rate of 0.25 (mean(c(0, 0.5)) rather than 0.5, as
before.