Solving for Juvenile Survival: Example Calculations with a Stage-Classified Population Matrix in R
John R. Brandon < jbrandon@gmail.com >
2015-08-14
Introduction
The goal of this document is to step through examples towards the calculation of a life history parameter value in a population dynamics model. The end product of this exercise will be the solution for the juvenile survival rate that would result in a pre-specifed population growth rate, holding constant other life history parameters (e.g., adult survival, birth rate, etc.). The approach taken here involves some principles from linear algebra and numerical methods for finding the root(s) of a function. For more background underlying an example motivation for this problem – including walrus and bowhead whales as case studies – see Brandon et al. (2007)(and references therein).
This document started as a working laboratory notebook. The text provides notes and snippets (aka “chunks”) of R code. It was initially drafted as a way to check solutions, by referencing results from R code, against a module of Fortran code that was being developed to preform the same task. I’ve since expanded the text in places, with a wider audience in mind.
This work is related to an open-source project involving a management strategy evaluation (see e.g. Punt et al. 2014). The ultimate goal of the project is to provide a platform for testing approaches that calculate allowable limits of human caused mortality, and meet the management objectives of the U.S. Marine Mammal Protection Act. For more background on the project, check out this webpage.
Step 1: Initialize a stage-classified matrix
A population dynamics model can be projected such that: \(\boldsymbol{n}(t + 1) = \boldsymbol{A} \boldsymbol{n}(t)\). Where: \(\boldsymbol{n}(t)\) is a vector of the numbers in each age/stage class at time t, and A is the projection matrix. The projection matrix is composed of birth and survival rates, and potentially other transitional rates in the life cycle (e.g. maturation or growth rates, etc…).
A stage-classified matrix “lumps” age-classes into demographic stages. This lumping collapses the size of the population projection matrix. If you are modeling the population dynamics of long-lived species, this can be a handy approach.
In this example, animals past a certain age (i.e. sexually mature individuals) will be lumped into a stage denoted by age x. The assumption being that vital rates are constant after a certain age (i.e. ignoring senescence). The model we will consider is actually a mix of a stage and ages. For brevity, I’ll refer to this as a stage-classified model (rather than, say, an age/stage-classified model). See the figure below to get a better sense of how age x is modeled (also known as a “plus-group” or “self-loop group”).
First we’ll initialize some life history parameters for a hypothetical population:
s_juv = 0.75 # Juvenile survival
s_adult = 0.85 # Adult (mature) survival
a_m = 4 # Age at sexual maturity
age_x = a_m + 1 # The "plus-group" stage (aka self-loop group)
b_max = 1.0 # Maximum annual birth rateFor more complicated applications we might instead read these values from an input file (e.g. using read.csv), or generate random variates from prior distributions if using Bayesian methods. For now though, we’ll limit the scope of this step and move on.
We’ll write a short function to assign values to the elements of the projection matrix A. This will make this step more easily repeatable as we continue to develope the code for this exercise in subsequent steps:
init_matrix = function(s_j, s_a, a_mat, age_xx, b_maxx){
A_matrix = matrix(data = 0, nrow = age_xx, ncol = age_xx) # dimension matrix and fill with zeros
A_matrix[1, (a_mat:age_xx)] = b_maxx # assign birth rates to mature ages in first row of matrix
for(ii in 1:(a_mat-1)) A_matrix[ii+1, ii] = s_j # assign juvenile survival rates
A_matrix[age_xx,a_mat] = s_a # adult survival assumed for maturing animals transitioning into plus-group
A_matrix[age_xx, age_xx] = s_a # adult survival assumed for plus-group
return(A_matrix)
}We can now call the function to assign values to the elements of the projection matrix, and print the matrix to check that our function is working correctly.
A = init_matrix(s_j = s_juv, s_a = s_adult, a_mat = a_m,
age_xx = age_x, b_maxx = b_max) # Call function to initialize matrix A
A # show matrix## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.00 0.00 0.00 1.00 1.00
## [2,] 0.75 0.00 0.00 0.00 0.00
## [3,] 0.00 0.75 0.00 0.00 0.00
## [4,] 0.00 0.00 0.75 0.00 0.00
## [5,] 0.00 0.00 0.00 0.85 0.85The figure below shows the life cycle graph that corresponds with the projection matrix A. Note that young of the year are often denoted as age zero in the literature, although R indexes arrays starting at one. This leads to a little book-keeping, so it’s something to keep in mind. The figure below purposely indexes starting at age zero, as a reminder. Question: How could this example life cycle be further simplified?
Step 2: Solve for the intrinsic population growth rate
We can solve for the intrinsic growth rate of this population (\(\lambda\)) by calculating the dominant real eigenvalue of the projection matrix A (see section 4.4. in Caswell (2001)). In general, calculating the eigenvalues of non-symetrical matrices is not a trivial task. Thankfully, the code to do this in R is trivial, and this is where we can take a moment to express our appreciation to the good people of the R Core Team for making the world a better place (R Core Team, 2013).
As a quick aside, R calls numerical procedures from the Linear Algebra PACKages (LAPACK) through wrapper functions like eigen() (which calls DGEEV(), a compiled Fortran procedure). In general, interfacing with the LAPACK procedures is much easier to do in languages like R, Python (e.g. using Numpy), MatLab, etc… than it is to do using a compiled language like C or Fortran. Rejoice!
Stubben and Millagan (2007) have provided the R package popbio which includes several helpful wrapper functions that simplify the (already relatively simple) R code necessary for interfacing with the LAPACK procedures involved in these calculations.
library("popbio") # Helpful wrapper functions (Stubben and Millagan, 2007)## Loading required package: quadproglambda(A) # Returns the dominant (largest in absolute value) real eigenvalue of the matrix, which is the intrinsic growth rate of the pop## [1] 1.137008So, for the life history parameter values we’ve assigned to this particular population dynamics model, the intrinsic rate of population growth would be 13.7% per year.
To see that the popbio function lambda() is returning the dominant (maximum absolute) real eigenvalue, we can type the function name (sans arguments) and inspect the code:
lambda## function (A)
## {
## ev <- eigen(A)
## lmax <- which.max(Re(ev$values))
## lambda <- Re(ev$values[lmax])
## lambda
## }
## <environment: namespace:popbio>So, lambda() is like a bow around the wrapper function eigen().
As another aside, the stable stage structure of the population can be calculated from the right eigenvector of the transition matrix A (see section 4.5 in Caswell (2001))
stable_age = stable.stage(A) #
stable_age## [1] 0.30943079 0.20410863 0.13463539 0.08880902 0.26301617Next let’s sum over the vector returned by stable.stage().
sum(stable_age)## [1] 1So, the function stable.age() returns the percentage of the population in each stage-class. For this example, 30.9% of the population will be young of the year in a stable stage distribution.
Step 3: Solve for the juvenile survival rate resulting in a specified value for \(\lambda\)
We want to solve for the juvenile survival rate that would result in a given intrinsic rate of population growth. For example, if we hold other life history parameters constant, what is the juvenile survival rate that would result in an intrinsic rate of population growth equal to 4% per year (\(\lambda = 1.04\))?
There are two approaches to solving this problem. The interested reader should note that there generally exists an analytical solution involving the characteristic equation of the matrix A, i.e. \(det(\boldsymbol{A} - \lambda \boldsymbol{I}) = 0\). Where: det is the detmerminant of the matrix; \(\lambda\) is a vector of eigenvalues (the roots of the characteristic equation); and I is the identity matrix.
Here, we’ll explore an alternative approach and tackle the problem using our computer and numerical algorithms. One advantage of numerical methods is that they can be easily re-used if the nature of the projection matrix changes. For example, if we wanted to model the population dynamics of plants, or invertebrates, we might use a size based model. A sized based model would probably not have the same characteristic equation as our relatively simple model here. So we’d have to go back to the drawing board and derive another analytical solution (and the det() can get messy quickly). With numerical methods, however, all we need to do is have a procedure that assigns values to the projection matrix and pass that info to the established coded procedures below. The remaining number crunching is generic. That’s nice.
In this final step, we’ll invoke a numerical minimization routine that will iteratively calculate values for \(\lambda\) over a range of juvenile survival rates, until it finds the solution for juvenile survival resulting in a target \(\lambda\).
First, we’ll need to define a function that calculates a value of \(\lambda\) resulting from a proposed value for juvenile survival. This \(\lambda\) value will be compared with the target value of 1.04. We’ll also need to define an objective function to be minimized. A natural choice for the objective function is the squared difference between the proposed and target value for \(\lambda\).
target_value = function(s_juv_tmp, s_a, a_mat, age_xx, b_maxx, target_lambda, dat.frame = TRUE){
A_matrix = init_matrix(s_juv_tmp, s_a, a_mat, age_xx, b_maxx)
lambda_tmp = lambda(A_matrix) # intrinsic rate of growth given value for s_juv_tmp
obj_fun = (lambda_tmp - target_lambda)^2 # objective function to be minimized
if(dat.frame){
return(data.frame(s_juv_tmp, lambda_tmp, obj_fun))
}else{
return(obj_fun)
}
}As an intermediate step, let’s take a look at the tabled relationship between juvenile survival and \(\lambda\) (again, holding other life history parameters constant).
param_vals = NULL
s_juv_range = seq(from = 0.5, to = 0.95, by = 0.01)
for(ii in 1:length(s_juv_range)){
param_vals = rbind(param_vals, target_value(s_juv_tmp = s_juv_range[ii],
s_a = s_adult, a_mat = a_m, age_xx = age_x,
b_maxx = b_max, target_lambda = 1.04))
}
head(param_vals, n = 20L) # Print first "n" param_values## s_juv_tmp lambda_tmp obj_fun
## 1 0.50 0.9820004 3.363959e-03
## 2 0.51 0.9876780 2.737595e-03
## 3 0.52 0.9934204 2.169662e-03
## 4 0.53 0.9992241 1.662675e-03
## 5 0.54 1.0050858 1.219004e-03
## 6 0.55 1.0110022 8.408730e-04
## 7 0.56 1.0169703 5.303666e-04
## 8 0.57 1.0229872 2.894354e-04
## 9 0.58 1.0290501 1.199012e-04
## 10 0.59 1.0351562 2.346209e-05
## 11 0.60 1.0413032 1.698262e-06
## 12 0.61 1.0474885 5.607721e-05
## 13 0.62 1.0537098 1.879593e-04
## 14 0.63 1.0599651 3.986032e-04
## 15 0.64 1.0662521 6.891710e-04
## 16 0.65 1.0725689 1.060733e-03
## 17 0.66 1.0789137 1.514273e-03
## 18 0.67 1.0852846 2.050695e-03
## 19 0.68 1.0916800 2.670821e-03
## 20 0.69 1.0980982 3.375406e-03Now, let’s inspect a couple of plots to make sure our calculations are sensisible, before moving onto the numerical minimization. First we’ll plot the relationship between juvenile survival and \(\lambda\), with a dashed horizontal line showing the \(\lambda\) target value.
library("ggplot2") # Hadley Wickham's excellent ggplot package for plotting
ggplot(data = param_vals, aes(x = s_juv_tmp, y = lambda_tmp)) + geom_line() +
xlab("Juvenile Survival Rate") + ylab("Intrinsic Pop Growth Rate (Lambda)") +
geom_abline(intercept = 1.04, slope = 0, lty = 2)
The relationship makes sense: Higher juvenile survival rates result in higher population growth rates, all else being equal. The solution for juvenile survival in this plot is where the two lines cross. Now, let’s take a look at the objective function (the squared difference between the calculated and targeted growth rate).
ggplot(data = param_vals, aes(x = s_juv_tmp, y = obj_fun)) + geom_line() +
xlab("Juvenile Survival Rate") + ylab("Objective Function Value")
From this plot, we can visualize the solution for the juvenile survival rate by inspecting the minimum of the objective function. We will now calculate a precise answer using the optimize function.
solution = optimize(f = target_value, interval = c(0.0, 1.0),
s_a = s_adult, a_mat = a_m, age_xx = age_x,
b_maxx = b_max, target_lambda = 1.04, dat.frame = FALSE)
solution## $minimum
## [1] 0.5978852
##
## $objective
## [1] 6.186673e-15There we have it, for this combination of life history parameters in this stage-classified population dynamics model, a juvenile survival rate of 0.5978852 corresponds with an intrinsic rate of population growth of 4% per year.
As a quick check:
new_matrix = init_matrix(s_j = solution$minimum, s_adult, a_m, age_x, b_max)
lambda(new_matrix)## [1] 1.04==================================== Happy modeling! ====================================