Life history trait practical using demographic models
Rob Salguero-Gómez (rob.salguero@biology.ox.ac.uk)
11/10/2025
Introduction to structured population models
In your lectures in year 1 and 2, you have been introduced to matrix population models (MPMs, hereafter) and -more briefly- to integral projection models (IPMs, hereafter). One of the commonalities of MPMs and IPMs is that they both describe the dynamics of a given population in discrete (rather than continuous) time. Another commonality is that, to describe those dynamics, MPMs and IPMs explicitly take into account how key demographic processes, such as survival and reproduction (a.k.a. “vital rates”), change among individuals in the study population. These models examine how individual heterogeneity relates to the values of these vital rates by explicitly studying the role that one (or more) individual-level characteristics (a.k.a. “states”) has on the vital rates. In the case of MPMs, individuals in the population are typically classified by age (e.g. 0, 1, 2 years old), developmental stage (e.g. juvenile, adult, senescent), and/or size ranges (e.g. short, medium, tall). Importantly here, when a researcher uses an MPM to parameterise the dynamics of individuals based on size, they must make a somewhat arbitrary call regarding where the boundary between, say, “short” vs. “medium”, and “medium” vs. “adult”, when individuals are classified along discrete size classes.
In this practical, you will use a dataset containing 1,000 individuals of the mystical golden horn unicorn (Monokeros aurum).
As you may remember, you went twice (in summer of 2019, and summer 2020) on a ecological expedition to the magical island of Monocerotum, where this species is endemic. Because the island is pretty small, and there are not hide-out places, the detectability of all unicorn individuals in the island is perfect. This feature of the island greatly simplifies how you will model this species’ population dynamics. The perfect detectability is of course also due to the fact that researchers at Oxford have been going to that island for over a decade, and the shy unicorns have acclimated to human presence, so they no longer turn invisible when they see a human. The fact that it is an island, and that unicorns are only found in it -and no where else around the globe- means that this is a “closed population”. As such, equation 1 below can be reduced to equation 2, thus simplifying ways in which we will model the golden unicorn population.
N(t+1) = N(t) + B + I - D - E (Equation 1)Just to remind you of the BIDE equation: N is the population size at a given time (t); B is the number of births within the population between time t and t+1; I is the number of individuals immigrating into the population between time t and t+1; D is the number of deaths in the population between time t and t+1; and E is the number of individuals emigrating from the population between time t and t+1.
N(t+1) = N(t) + B - D (Equation 2)It is important to note that equations 1 and 2 describe the dynamics of a population where individual heterogeneity has not been taken into account. Furthermore, equation 2 assumes a closed population (i.e. no immigration, no emigration), as in the case of the examined golden unicorn population.
Question 1 (Q0). Briefly explain in biological terms equation 1.
The “collection” of data
As you may hopefully remember, you went to the island of Monocerotum in the summer of 2019 and summer of 2020. In 2019, you sampled all existing 1,000 individuals of the golden unicorn (fun fact: this species is asexual, and so there are no males and females). For each individual, you measured the length of the horn with a measuring tape. We will call this variable ‘horn size’ from here on. You also tagged the right ear of each unicorn in 2019. The tags allowed you to examine in 2020 which individuals survived between 2019 and 2020, as opposed to individuals that were born in that time interval, which of course do not yet have an ear tag. In the following code chunk, you will be able to explore and plot some basic relationships between the different variables that you measured in the island of Monocerotum for this mystical species.
#Read in the dataset. Note that you will have to first change the writing directory (with the command "setwd(...)") to wherever you saved the csv file
data <- read.csv("GoldenUnicornData.csv", nrows = FALSE)
#Reading the first 20 rows of the data
head(data, n = 20) ## ID Horn.t0 Survival Reproduction Horn.t1 OffspringHorn.t1
## 1 1 9.30 0 0 NA NA
## 2 2 11.45 0 0 NA NA
## 3 3 12.45 0 0 NA NA
## 4 4 15.95 0 0 NA NA
## 5 5 16.30 0 0 NA NA
## 6 6 16.45 0 0 NA NA
## 7 7 17.85 0 0 NA NA
## 8 8 18.15 0 0 NA NA
## 9 9 18.40 1 0 19.4 NA
## 10 10 18.95 0 0 NA NA
## 11 11 19.35 0 0 NA NA
## 12 12 20.00 0 0 NA NA
## 13 13 20.15 0 0 NA NA
## 14 14 20.55 0 0 NA NA
## 15 15 20.85 0 0 NA NA
## 16 16 20.90 0 0 NA NA
## 17 17 21.60 0 0 NA NA
## 18 18 21.60 1 0 20.6 NA
## 19 19 22.10 0 0 NA NA
## 20 20 22.45 0 0 NA NAThis dataset is organised by descending horn size value in 2019, from the smallest value (9.30 cm) to the largest value (see below, 72.55 cm)
#Reading the last 20 rows of the data
tail(data, n = 20) ## ID Horn.t0 Survival Reproduction Horn.t1 OffspringHorn.t1
## 981 981 62.80 0 1 NA 16.6
## 982 982 62.90 1 1 66.40 18.6
## 983 983 63.00 1 1 63.50 17.2
## 984 984 63.20 1 1 65.70 17.2
## 985 985 63.40 0 0 NA NA
## 986 986 63.45 1 1 64.95 20.2
## 987 987 63.70 1 1 67.20 20.2
## 988 988 63.90 1 0 66.90 NA
## 989 989 63.95 0 0 NA NA
## 990 990 64.30 1 1 67.30 20.2
## 991 991 64.50 1 0 64.50 NA
## 992 992 65.15 0 0 NA NA
## 993 993 65.85 1 0 67.85 NA
## 994 994 66.35 0 0 NA NA
## 995 995 67.60 1 0 70.60 NA
## 996 996 67.85 1 0 68.35 NA
## 997 997 68.55 0 0 NA NA
## 998 998 69.25 1 0 71.25 NA
## 999 999 71.30 0 0 NA NA
## 1000 1000 72.55 1 0 73.55 NAThe variables in this dataset are:
- ID: unique identifier for each unicorn, corresponding to the unique ear tag you attached in 2019
- Horn.t0: length of the horn in 2019, in cm
- Survival: whether the individual survived (1) or not (0) from 2019 to 2020
- Reproduction: whether the individual reproduced (1) or not (0) in 2019 and its offspring survived till you sampled it again in 2020
- Horn.t1: length of the horn in 2020, in cm
- OffspringHorn.t1: length of the horn of the new bundles of joy (a.k.a offspring) in 2020
The organisation of the data above show a few logical, biological aspects:
- Check out ID 9, 18, 982, 983, etc: only individuals that survive to 2020 have a horn size (i.e different than NA) in 2020.
- Check out ID 982, 986, 990: only individuals that reproduce in 2019 can produce an offspring of a given horn size in 2020
The information on survival (the opposite to death - “D” in equations 1, 2) and reproduction (“B” in equations 1, 2) allow you to describe the basic dynamics of the population. However, since individuals vary quite a bit in their vital rates as a function of their horn size, it may be worth first exploring whether horn size helps predict survival and reproduction. In other words: do individuals in the population show a great deal of individual heterogeneity in vital rates as a function of their horn size? To answer this question, let us explore the dataset with a few graphs:
- The relationship between horn size in 2019 and survival to 2020:
plot(data$Horn.t0, jitter(data$Survival, 0.2), xlim = c(0, 75), ylab = "Survival probability", xlab = "Horn length in 2019 (cm)") #the command "jitter", as used here, spreads the data a bit for visualisation purposes along the y axis, so you can see how individuals of a give horn size may vary in their survival outcomes between 2019 and 2020
Question 1 (Q1). Describe the relationship in the figure above in biological terms.
- The relationship between horn size in 2019 and reproduction in 2019:
plot(data$Horn.t0, jitter(data$Reproduction, 0.2), xlim = c(0, 75), ylab = "Reproduction probability", xlab = "Horn length in 2019 (cm)") #jitter is used here in a similar way as above for reproduction
Q2. Describe the relationship in the figure above in biological terms.
- The relationship between horn size in 2019 and the size of the horn of the same individual in 2020:
plot(data$Horn.t0, data$Horn.t1, xlim = c(0, 75), ylim = c(0, 75), ylab = "Horn length of surviving individual in 2020 (cm)", xlab = "Horn length in 2019 (cm)" )
abline(0, 1, lty = 2, col = "gray")
Q3. Describe the relationship in the figure above in biological terms. What does the dashed, gray line represent in biological terms?
Q4. Now, and because you are so interested in the biology of the golden unicorn, plot the relationship (and present the fully commented R script to do so) between the horn size of reproductive individuals in 2019 and horn size of their offspring in 2020. Does there seem to be a relationship? If so, what biological mechanisms could explain it?
#Make sure you also insert the correct labels on the x and y axes!
#Hint: Take inspiration from the chunks above and re-purposes commands to fit your needsDescribe the relationship in the figure above in biological terms.
From the figures above, you can se that horn size in 2019 should be (i.e. pending statistical tests, below) a good predictor of survival and reproduction. Thus, it would be sensible to build a model where horn size is a state predictor. As the data are discrete in time (e.g. you have data for 2019 and 2020, but not in continuous manner), we can chose matrix population models or integral projection models to model the dynamics of that population. The following section will show you how to build a MPM.
Building your Matrix Population Model
As you may remember from recent lectures, an MPM is a discrete time, discrete state demographic model. The data you have collected, however, have classified individuals along a continuous trait: horn size. Thus, the first key decision is to categorise individuals into discrete classes. Here, we will chose three classes and call them “A” (for individuals with horns shorter than 30 cm), “B” (for individuals with horn length between 30 and 60), and “C” (for individuals with horns longer than 60 cm), as described by:
- A: < 30
- B: >30 & < 60
- C: >60
In the following chunk of code, we will assign individuals to each of these three stages in 2019 and 2020:
#Create new variables where the data will be saved to discretise classes. For the time being we will assign them the empty value "NA"
data$Horn.Discrete.t0 <- NA
data$Horn.Discrete.t1 <- NA
data$OffspringHorn.Discrete.t1 <- NA
#Define the thresholds that separate class A from B:
smallBoundary <- 30
#Define the thresholds that separate class B from C:
largeBoundary <-60
#Assign sizes to discrete classes for individuals measured in 2019:
data$Horn.Discrete.t0[which(data$Horn.t0 < smallBoundary)] <- "A" #We start by assigning individual with horns smaller than 30 cm to class "A"
data$Horn.Discrete.t0[which(data$Horn.t0 >= smallBoundary & data$Horn.t0 <= largeBoundary)] <- "B" #Next we assign individuals with horns larger than 30 and smaller than 60 to class "B"
data$Horn.Discrete.t0[which(data$Horn.t0 > largeBoundary)] <- "C" # finally we assign individuals with horns larger than 60 to class "C"
#Let's now plot the raw data for 2019 and take a look at where you placed the thresholds
hist(data$Horn.t0)
abline(v = smallBoundary, col = "purple", lwd = 3, lty = 2) #Adds a "v"ertical line to the plot at a given value, in this case smallBoundary
#Now insert a vertical line showing in "orange" colour the largeBoundary value in the plot above
#hint: abline(v = ...)
#Inspect how many individuals now fall in each of your classes in 2019:
table(data$Horn.Discrete.t0)##
## A B C
## 95 865 40Q5. Describe the histogram above.
#Do the same thing as above for size in 2020
data$Horn.Discrete.t1[which(data$Horn.t1 < smallBoundary)] <- "A"
data$Horn.Discrete.t1[which(data$Horn.t1 >= smallBoundary & data$Horn.t1 <= largeBoundary)] <- "B"
data$Horn.Discrete.t1[which(data$Horn.t1 > largeBoundary)] <- "C"
#And plot it again for 2019
hist(data$Horn.t1)
abline(v = smallBoundary, col = "purple", lwd = 3, lty = 2)
abline(v = largeBoundary, col = "orange", lwd = 3, lty = 2)
#Inspect how many individuals now fall in each of your classes in 2020:
table(data$Horn.Discrete.t1)##
## A B C
## 17 730 40Q6. Describe the histogram above.
#Finally we classify offspring too into these way of discretising size
data$OffspringHorn.Discrete.t1[which(data$OffspringHorn.t1 < smallBoundary)] <- "A"
data$OffspringHorn.Discrete.t1[which(data$OffspringHorn.t1 >= smallBoundary & data$OffspringHorn.t1 <= largeBoundary)] <- "B"
data$OffspringHorn.Discrete.t1[which(data$OffspringHorn.t1 > largeBoundary)] <- "C"
#And plot it again for 2020
hist(data$OffspringHorn.t1)
abline(v = smallBoundary, col = "purple", lwd = 3, lty = 2)
abline(v = largeBoundary, col = "orange", lwd = 3, lty = 2)
#Inspect how many offspring fall in each of your classes in 2020:
table(data$OffspringHorn.Discrete.t1)##
## A
## 281Q7. Describe the histogram above. Why do your vertical lines not show?
Now that your data are discretised, we can compile the MPM using a series of standard steps:
Step 1. Quantify the proportion of individuals from each of the discretised size classes that survived and transitioned to each of the other three classes
#This matrix shows how many surviving individuals in each class transition to each class in the next year:
matrixTransitions <- table(data$Horn.Discrete.t1,data$Horn.Discrete.t0)
matrixTransitions <- matrix(data= matrixTransitions, nrow = length(unique(data$Horn.Discrete.t0)), ncol= length(unique(data$Horn.Discrete.t0))) #Converting this table into a matrix so that calculations we will do with R below will work well.
matrixTransitions## [,1] [,2] [,3]
## [1,] 17 0 0
## [2,] 11 719 0
## [3,] 0 12 28Q8. Describe what the numbers in the matrix above mean in biological terms.
Step 2. Now convert the frequency matrix to probabilities of transition:
#Calculate how many individuals are transitioning from the different classes in 2019
classTransitions <- colSums(matrixTransitions) #The command colSums sums the columns of a matrix
classTransitions## [1] 28 731 28#Convert the frequency of each column in the matrix into a probability of transitioning to different classes in 2020
matrixTransitions[,1] <- matrixTransitions[,1]/classTransitions[1]
matrixTransitions[,2] <- matrixTransitions[,2]/classTransitions[2]
matrixTransitions[,3] <- matrixTransitions[,3]/classTransitions[3]
matrixTransitions## [,1] [,2] [,3]
## [1,] 0.6071429 0.00000000 0
## [2,] 0.3928571 0.98358413 0
## [3,] 0.0000000 0.01641587 1Q9. Describe what the numbers in the matrix above mean in biological terms.
Of course, this matrix does not take into account the individuals that were alive in 2019, but did not make it to 2020, so it needs to be “penalised” by the stage-specific survival probabilities. For that:
Step 3. Penalise the different transition probabilities by class-specific survival values
#This matrix shows how many individuals in each class transition to each class in the next year. Remember: 0 means death, 1 means survival
classSurvival <- table(data$Survival,data$Horn.Discrete.t0)
classSurvival##
## A B C
## 0 67 134 12
## 1 28 731 28#This matrix shows how many individuals in each class transition to each class in the next year:
classSurvival <- classSurvival[2,]/colSums(classSurvival)
classSurvival## A B C
## 0.2947368 0.8450867 0.7000000Q10. Describe what the numbers in the matrix above mean in biological terms.
#This matrix is now penalised by survival
matrixTransitions[,1] <- matrixTransitions[,1]*classSurvival[1]
matrixTransitions[,2] <- matrixTransitions[,2]*classSurvival[2]
matrixTransitions[,3] <- matrixTransitions[,3]*classSurvival[3]
matrixTransitions## [,1] [,2] [,3]
## [1,] 0.1789474 0.00000000 0.0
## [2,] 0.1157895 0.83121387 0.0
## [3,] 0.0000000 0.01387283 0.7Step 4. Construct a matrix that links the number of offspring that fall in each of the discrete classes to per-capita events of reproductive unicorns in 2019, as a function of their class
totalOffspringPerClass <- table(data$OffspringHorn.Discrete.t1, data$Horn.Discrete.t0)
totalOffspringPerClass##
## A B C
## A 0 260 21Q11. Describe what the numbers in the table above mean in biological terms.
Step 5. It is important to note that every time a unicorn gives birth, its gives birth to just one new cute little bundle of joy (and no more). Thus, all we need to do is to calculate the probability that individuals in each class reproduced is a table that tells us how many individuals reproduced (or not) as a function of their class in 2019
classReproduction <- table(data$Reproduction, data$Horn.Discrete.t0)
classReproduction##
## A B C
## 0 95 604 19
## 1 0 260 21And then we need to calculate the per-capita contributions as it follows:
classReproduction <- classReproduction[2,]/colSums(classReproduction)
classReproduction## A B C
## 0.0000000 0.3009259 0.5250000Q12. Describe what the numbers above mean in biological terms.
Step 6. As we know that those reproductive events go to class “A” (given that this class was originally defined to contain the largest of all new offspring too - see histograms above), we insert that vector in the first row of a reproduction matrix that has the same dimensions (i.e. same number of rows and columns) as the transition matrix containing the transition probabilities:
matrixReproduction <- matrix(0, length(unique(data$Horn.Discrete.t0)), length(unique(data$Horn.Discrete.t0)))
matrixReproduction[1,] <- classReproduction
MPM <- matrixTransitions + matrixReproduction
MPM## [,1] [,2] [,3]
## [1,] 0.1789474 0.30092593 0.525
## [2,] 0.1157895 0.83121387 0.000
## [3,] 0.0000000 0.01387283 0.700Q13. Which of the matrix elements above correspond to transition probabilities?
Q14. Which of the matrix elements above correspond to per-capita reproductive contributions?
Analysing your Matrix Population Model
Now that we have the model created, let us obtain some key demographic properties from it:
- The rate of growth of the population is:
library(popbio) #You will need to install this package if ou haven't already using the command "install.packages("popbio")"
lambda <- lambda(MPM)
lambda## [1] 0.886815Q15. Explain the value of lambda in biological terms.
- The generation time is:
library('Rage') #To install this package, follow the "Installation" instructions at the bottom of this link: https://github.com/jonesor/Rage
generationTime <- gen_time(matU = matrixTransitions, matR = matrixReproduction)
generationTime## [1] 10.84755Q16. Explain the value of generation time in biological terms.
- The elasticities of growth of the population with respect to the different matrix elements are:
library(popbio)
elasticities <- elasticity(MPM)
elasticities## [,1] [,2] [,3]
## [1,] 0.01424467 0.049885209 0.006462861
## [2,] 0.05634807 0.842379851 0.000000000
## [3,] 0.00000000 0.006462861 0.024216481Q17. Explain what an elasticity is in the context of an MPM.
#image(MPM)
library('plot.matrix') #install this package too if you don't have it already in your computer
par(mar = c(2, 2, 1, 3))
par(mgp = c(3, 1, 0))
par(oma = c(2, 1, 1, 1))
plot(elasticities)
Q18. Identify in the MPM above the matrix element with highest elasticity - explain in biological terms what that element represents.
- Mean life expectancy
library('Rage') #To install this package, follow the "Installation" instructions at the bottom of this link: https://github.com/jonesor/Rage
MeanLifeExpectancy <- longevity(matU = matrixTransitions)
MeanLifeExpectancy## [1] 17Q19. Interpret in biological terms the value obtained above for mean life expectancy. What time units is this in? Why does the function ‘longevity’ not need the argument ‘matR = matrixReproduction’, as it was required for the calculatin of generation time above?
Q20. Up to here you have constructed and analysed a MPM based on continuous trait data (e.g. horn size) that was discretised following a series of arbitrary criteria/thresholds that help separate indidivuals into three classes. A question that (I hope!) would have already popped in your brain is: ‘what happens if I change those arbitrary thresholds? How sensitive are the outputs to my decisions?’. Next, please answer Questions 5-19 under the following values of thresholds for the different classes: - Class A: individuals with horn size < 25 - Class B: between 25 and 50 - Class C: > 50 Did the results in sections a and b for Q 5-19 change. If so, why do you think that is?
Luckily, there are models that allow you to quantify the dynamics of a population while taking into consideration individual heterogeneity based on continuous trait data (e.g. as here with horn size). These models are called integral projection models (IPMs) and, while they’ve been briefly introduced in lecture, they offer a series of advantages over MPMs that may be of your interest going forward. Some of these include them being more robust to small sample sizes, in addition to them not needed pre-fixed, arbitrary thresholds to defined discrete classes. Below, I have include a few references for MPM and IPM modelling.
Further references
The following citations are of interest for this practical:
https://www.jstor.org/stable/177363?seq=1
https://www.sinauer.com/media/wysiwyg/tocs/MatrixPopulationModels2.pdf
In addition, I also note that, because the horn data in this dataset is in fact a continuous trait, these dataset is highly amenable to the construction of integral projection models (IPMs). I am happy to run tutorials on IPMs on request in Hillary/Trinity. In the meantime, some important references on IPMs include:
https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/0012-9658(2000)081[0694:SSSAAN]2.0.CO;2
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210x.12001
https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.12178
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12146
https://www.springer.com/gp/book/9783319288918
Please note that the last resource is accessible open-access online/via the university library system.