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.

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

You can refer to the documentation for the lifeTimeRepEvents function in the Rage package to better understand all of the arguments, inputs, and outputs of this function here.

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
#> [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?