Last week in Davis R Users’ Group, Lauren Yamane showed us how she created and analyzed a stochastic age-structured population in R. Her examples are below. Her original scripts can be found as *.Rmd files here

A note to UC Davis students: This topic and others will be covered by Marissa Baskett and Sebastian Schreiber in their course this winter, Computational methods in population biology (ECL298)


A discrete time, age-structured model of a salmon population (semelparous) that can live to age 5, with fishing and environmental stochasticity

Parameter values

Make a function for the age-structured matrix

Run model by calling function

Plot of time series with all 5 age classes

Export text file for plot with age classes

Or plot of time series for total population

Export text file for plot with total population

Stability analysis of deterministic population model

This assumes that the partial differential equations that make up the Jacobian matrix have been determined analytically. Below are the parameters that were associated with my age-structured salmon model with fishing. I can compare the stability of the population with fishing to that with no fishing

Fishing Jacobian

#         [,1] [,2]   [,3]     [,4]  [,5]
# one_F   0.00 0.00 2.5214 21.35558 23.88
# two_F   0.28 0.00 0.0000  0.00000  0.00
# three_F 0.00 0.28 0.0000  0.00000  0.00
# four_F  0.00 0.00 0.2504  0.00000  0.00
# five_F  0.00 0.00 0.0000  0.02957  0.00

Fishing Eigenvalues

No Fishing Jacobian

#           [,1] [,2]     [,3]    [,4]   [,5]
# one_noF   0.00 0.00 0.004625 0.03917 0.0438
# two_noF   0.85 0.00 0.000000 0.00000 0.0000
# three_noF 0.00 0.85 0.000000 0.00000 0.0000
# four_noF  0.00 0.00 0.760240 0.00000 0.0000
# five_noF  0.00 0.00 0.000000 0.08976 0.0000

No Fishing Eigenvalues

Plot the eigenvalues with legend to compare stability of both deterministic models

layout(matrix(c(1, 0, 2, 0, 3, 0), 2, 3, byrow = TRUE))  #Layout with 2 rows, 3 columns,
# the matrix tells you the number of the figure in each location. Since
# the layout is byrow, it fills the first row first and then the 2nd row
# no fishing plot
plot(ev1_noF[1], ev1_noF[2], type = "p", pch = 0, cex = 2, xlab = "Real", ylab = "Imaginary", 
    xlim = c(-1.25, 1.25), ylim = c(-1.25, 1.25))
title("No fishing", adj = 0)  #plot dominant eigenvalue, with x-coordinate as real part and y-coordinate as imaginary part
draw.circle(0, 0, radius = 1, border = "black", lty = 1, lwd = 1)  #draw unit circle
points(ev2_noF[1], ev2_noF[2], pch = 1, cex = 2)  #2nd eigenvalue
points(ev3_noF[1], ev3_noF[2], pch = 2, cex = 2)  #3rd eigenvalue
points(ev4_noF[1], ev4_noF[2], pch = 6, cex = 2)  #4th eigenvalue
points(ev5_noF[1], ev5_noF[2], pch = 5, cex = 2)  #5th eigenvalue
abline(a = NULL, b = NULL, h = 0, lty = 3)
abline(a = NULL, b = NULL, v = 0, lty = 3)

# fishing plot
plot(ev1_F[1], ev1_F[2], type = "p", pch = 0, cex = 2, xlab = "Real", ylab = "Imaginary", 
    xlim = c(-1.25, 1.25), ylim = c(-1.25, 1.25))
title("Fishing", adj = 0)
draw.circle(0, 0, radius = 1, border = "black", lty = 1, lwd = 1)
points(ev2_F[1], ev2_F[2], pch = 1, cex = 2)
points(ev3_F[1], ev3_F[2], pch = 2, cex = 2)
points(ev4_F[1], ev4_F[2], pch = 6, cex = 2)
points(ev5_F[1], ev5_F[2], pch = 5, cex = 2)
abline(a = NULL, b = NULL, h = 0, lty = 3)
abline(a = NULL, b = NULL, v = 0, lty = 3)

# Legend
plot(1:10, type = "n", axes = FALSE, xlab = "", ylab = "", frame.plot = TRUE)
title("Legend", adj = 0)
text(6.5, 9.5, expression(lambda[1]), cex = 1.2)
text(6.5, 7.5, expression(lambda[2]), cex = 1.2)
text(6.5, 5.5, expression(lambda[3]), cex = 1.2)
text(6.5, 3.5, expression(lambda[4]), cex = 1.2)
text(6.5, 1.5, expression(lambda[5]), cex = 1.2)
points(4.5, 9.5, pch = 0, cex = 1.3)
points(4.5, 7.5, pch = 1, cex = 1.3)
points(4.5, 5.5, pch = 2, cex = 1.3)
points(4.5, 3.5, pch = 6, cex = 1.3)
points(4.5, 1.5, pch = 5, cex = 1.3)


← A multi-species, age-structured model for forest disease | All posts | Qualifying Exam Reading List →