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.

Walkthrough and R code

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)
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) {
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)-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

1. Individuals can grow, shrink, and stay the same size
2. Individuals can reproduce sexually
3. 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
#>  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
#>  3.974702

dyn\$stable.stage
#>  0.845445401 0.071542444 0.073388842 0.009623313

dyn\$repro.value
#>   1.000000  1.511156 11.404517  8.658467

N <- solve(diag(dim(matU))-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)
#>  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
#>  74.70513

generation.time
#>  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?