a02_VitalRates.Rmd

---
title: "Deriving vital rates from an MPM"
author: "Patrick Barks"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Deriving vital rates from an MPM}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```



## Introduction

The 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).

## Preliminaries: Constructing an MPM from lower-level vital rates

```{r}
library(Rage)
```

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.

```{r}
## 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 # large
```

Next we'll use the vital rates to construct an MPM, based on a post-breeding census design.

```{r}
# growth/survival component
matU <- rbind(
  c(s1 * g11, 0, 0, 0),
  c(s1 * g12, s2 * g22, 0, 0),
  c(s1 * g13, s2 * g23, s3 * g33, s4 * g43),
  c(0, s2 * g24, s3 * g34, s4 * g44)
)

# sexual reproduction component
matF <- rbind(
  c(0, s2 * f2, s3 * f3, 0),
  c(0, 0, 0, 0),
  c(0, 0, 0, 0),
  c(0, 0, 0, 0)
)
```

## Deriving vital rates of survival

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()`.

```{r}
(surv <- vr_vec_survival(matU))
# equivalent to...
(surv <- colSums(matU))
```

As expected, these values match the survival probabilities that we initially specified.

## Deriving vital rates conditional on survival

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.

```{r}
# matU vital rates conditional on survival
vr_mat_U(matU)

# matF vital rates conditional on survival
vr_mat_R(matU, matF)
```

We've now recovered all the vital rates that we initially specified when we constructed our MPM.

## Summarizing vital rates within stages

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.

#### Growth, shrinkage, stasis, and dormancy

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`.

```{r}
vr_vec_growth(matU)
```

However, 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.

```{r}
vr_vec_growth(matU, exclude_row = 4)
```

Next 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.

```{r}
vr_vec_shrinkage(matU)
```

However, 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.

```{r}
vr_vec_shrinkage(matU, exclude_col = 4)
```

Next 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.

```{r}
vr_vec_stasis(matU)
```

Finally, 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.

```{r}
vr_vec_dorm_enter(matU, dorm_stages = 4)
vr_vec_dorm_exit(matU, dorm_stages = 4)
```

#### Fecundity

The 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).

```{r}
vr_vec_reproduction(matU, matF)
```

In 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.

## Summarizing vital rates across stages

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.

#### Simple average across stage classes

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.

```{r}
vr_growth(matU, exclude_row = 4)

# equivalent to
(vec_growth <- vr_vec_growth(matU, exclude_row = 4))
mean(vec_growth, na.rm = TRUE)
```

Note 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.

#### Weighted average across stage classes

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.

```{r}
# calculate the stable distribution using popdemo::eigs
library(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)

# equivalent to
pos <- !is.na(vec_growth) # possible transitions
sum(vec_growth[pos] * w[pos]) / sum(w[pos]) # weighted average
```

## Distinguishing between possible and impossible transitions

By 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).

```{r}
vr_vec_shrinkage(matU)
vr_shrinkage(matU)
```

However, 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.

```{r}
(pos <- matU > 0) # possible transitions based on matU
pos[2, 3] <- TRUE # set medium-to-small shrinkage tr to possible

vr_vec_shrinkage(matU, posU = pos)
vr_shrinkage(matU, posU = pos)
```

Now, 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.

## References

Caswell, H. (2001). Matrix Population Models: Construction, Analysis, and Interpretation. 2nd edition. Sinauer Associates, Sunderland, MA. ISBN-10: 0878930965