Introduction

In addition to population growth rate and stable stage distributions, MPMs can be used to calculate a number of other important life history metrics. Many of these calculations are included as functions in the *Rage* package, but here we’ll show you the underlying code for the calculations in R.

Video tutorial

A video tutorial has been created to accompany this page and demonstrate these excercises.

Check it out here on the COM(P)ADRE Learn YouTube channel!

Walkthrough and R code

The R code for these exercises can be downloaded as an R file or PDF in addition to being displayed on this page.

Download as R file

Download as PDF

Preliminaries

For this exercise, we’ll use the *MASS* and *popbio* packages. We’ll reference the *Rage* package, which is in progress and will likely be released before summer 2021.

```
``````
library(MASS)
library(popbio)
```

Other key life history traits from MPMs

We will start by writing a function to calculate life history traits from a matrix. Specifically, the traits we will examine are the probability of surviving to first reproductive event, age at maturity, and reproductive lifespan. This code is primarily adapted from **matlab** code from Hal Caswell’s textbook on matrix population models (Caswell 2001).

```
``````
lifeTimeRepEvents <- function(matU, matF, startLife = 1) {
uDim = dim(matU)[1]
surv = colSums(matU)
repLifeStages = colSums(matF)
repLifeStages[which(repLifeStages>0)] <- 1
if(missing(matF) | missing(matU)) {stop('matU or matF missing')}
if(sum(matF, na.rm=T)==0) {stop('matF contains only 0 values')}
# Probability of survival to first reprod event
Uprime <- matU
Uprime[, which(repLifeStages==1)] <- 0
Mprime = matrix(0,2,uDim)
for (p in 1:uDim[1]) {
if (repLifeStages[p]==1) {Mprime[2,p] = 1} else {
Mprime[1,p] = 1-surv[p]
}
}
Bprime = Mprime%*%(ginv(diag(uDim)-Uprime))
pRep = Bprime[2,startLife]
out = data.frame(pRep = pRep)
# Age at first reproduction (La; Caswell 2001, p 124)
D = diag(c(Bprime[2,]))
Uprimecond = D%*%Uprime%*%ginv(D)
expTimeReprod = colSums(ginv(diag(uDim)-Uprimecond))
La = expTimeReprod[startLife]
out$La = La
# Mean life expectancy conditional on
# entering the life cycle in the first reproductive stage
firstRepLifeStage = min(which(repLifeStages==1))
N = solve(diag(uDim[1])-matU)
meanRepLifeExpectancy = colSums(N)[firstRepLifeStage]
out$meanRepLifeExpectancy = meanRepLifeExpectancy
# Life expectancy from mean maturity
remainingMatureLifeExpectancy = colSums(N)[startLife]-La
out$remainingMatureLifeExpectancy = remainingMatureLifeExpectancy
return(out)
}
```

At this point, we can either choose a matrix from COMPADRE/COMADRE or create our own matrix model. Here, we create a matrix where

- Individuals can grow, shrink, and stay the same size
- Individuals can reproduce sexually
- The population has a seedbank as its first stage

The matrix model will comprise two component matrices. The **U** matrix contains the survival-dependent processes or transitions. We will check the stage-specific survival rates, which cannot be > 1, by summing the columns in the **U** matrix. The **F** matrix contains the sexual reproduction rates, and the full (**A**) matrix is the sum of the **U** and **F** matrices (assuming no clonal reproduction for this population).

```
``````
matU <- matrix(c(0.1,0.1,0.3,0.1,0.3,0.3,0.1,0.2,0.3,0.3,0.2,0.2,0,0.2,0.3,0.2),
nrow = 4, ncol = 4, byrow = T)
colSums(matU) # No column can be > 1
#> [1] 0.7 0.9 0.9 0.7
matF <- matrix(c(0,0.3,40,30, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
nrow = 4, ncol = 4, byrow = T)
matA <- matU + matF
```

Let’s start calculating life history traits using the matrix model. First, we’ll take a quick look at the overall population dynamics based on the matrix using the `eigen.analysis`

function from the *popbio* package. Then, we’ll calculate the fundamental matrix, **N**, which gives the residence time of individuals in each of the stages and transitions. Summing each column in matrix **N** gives the mean life expectancy conditional on entering the life cycle in each of the stages.

```
``````
dyn <- eigen.analysis(matA)
dyn$lambda1
#> [1] 3.974702
dyn$stable.stage
#> [1] 0.845445401 0.071542444 0.073388842 0.009623313
dyn$repro.value
#> [1] 1.000000 1.511156 11.404517 8.658467
N <- solve(diag(dim(matU)[1])-matU)
N
#> [,1] [,2] [,3] [,4]
#> [1,] 1.7521368 0.892094 1.030983 0.6997863
#> [2,] 1.1538462 2.355769 1.105769 1.0096154
#> [3,] 1.2820513 1.506410 2.339744 1.1217949
#> [4,] 0.7692308 1.153846 1.153846 1.9230769
colSums(N)
#> [1] 4.957265 5.908120 5.630342 4.754274
```

Generation time (**T**) is calculated using the net reproductive rate (**Ro**) and lambda.

```
``````
R <- matF %*% N
Ro <- lambda(R)
lambda <- lambda(matA)
generation.time <- log(Ro)/log(lambda)
Ro
#> [1] 74.70513
generation.time
#> [1] 3.125874
```

Finally, we’ll use the function that we created above to estimate other life history traits. Here, we specify the option `startLife = 2`

because in our matrix model the first stage is for the seedbank, and typically in plant demographics, life starts after germination.

```
``````
lifeTimeRepEvents(matU, matF, startLife = 2)
#> pRep La meanRepLifeExpectancy remainingMatureLifeExpectancy
#> 1 1 1 5.90812 4.90812
```

**Extra:** Try changing this option to `startLife = 1`

to see what happens. Why do the results change, and what do the changes mean biologically?