This document walks through an elasticity analysis for a songbird population. Elasticity analysis is a matrix population modeling tool that tells us which vital rates have the greatest proportional influence on population growth rate. In applied contexts this information may inform management priorities, but the analysis itself describes ecological leverage within the life cycle.
This tutorial follows a straightforward workflow:
The three focal vital rates in this analysis are:
A matrix population model (also called a Leslie matrix or stage-based matrix model) is a mathematical framework that projects how many individuals will be in each age or stage class in the next time step, given the current number of individuals at each stage.
In this model, we track two stage classes:
We use a theoretical post-breeding census, meaning we count individuals after the breeding season each year. This is standard for many songbird studies because it lets us include new fledglings in the count immediately.
The model advances in annual time steps, so every quantity below represents one full year.
The projection matrix \(\mathbf{A}\) is a 2×2 table where each cell describes a demographic transition or contribution between stage classes:
\[ \mathbf{A} = \begin{pmatrix} 0 & F \\ \phi_j & \phi_a \end{pmatrix} \]
Reading this table:
| Cell | Symbol | Meaning |
|---|---|---|
| [1, 1] | 0 | Juveniles cannot yet reproduce — they contribute no new fledglings |
| [1, 2] | \(F\) | Fecundity — the per-capita contribution of adults to new juveniles |
| [2, 1] | \(\phi_j\) | Post-fledgling survival — probability a juvenile survives to become an adult |
| [2, 2] | \(\phi_a\) | Adult survival — probability an adult survives to breed again |
Columns represent the current stage and rows represent the future stage, so each column describes what happens to individuals in that stage over one year. Here is how we can visualize this matrix:
Note: There is no [1, 1] in this model because juveniles cannot reproduce!
The matrix is just a compact mathematical way of describing who survives and who produces new young.
# Point estimates — replace with species-specific values:
v_rate <-
list(
breeding_attempts = 1.8,
clutch_size = 3.5,
sex_ratio = 0.5,
nest_success = 0.35,
phi_j = 0.25,
phi_a = 0.55
)| Symbol | Parameter | Value |
|---|---|---|
| $$B$$ | Breeding attempts per season | 1.80 |
| $$CS$$ | Mean clutch size | 3.50 |
| $$SR$$ | Sex ratio | 0.50 |
| $$NS$$ | Nest success | 0.35 |
| $$\phi_j$$ | Post-fledgling annual survival | 0.25 |
| $$\phi_a$$ | Adult annual survival | 0.55 |
Fecundity \(F\) is the product of several biological rates that together determine how many new fledglings each adult female contributes to the population per year:
\[ F = B \times CS \times SR \times NS \]
where:
| Symbol | Parameter | Default value |
|---|---|---|
| \(B\) | Breeding attempts per season | 1.8 |
| \(CS\) | Mean clutch size (eggs per attempt) | 3.5 |
| \(SR\) | Sex ratio — proportion of offspring that are female | 0.5 |
| \(NS\) | Nest success — probability a nest produces ≥ 1 fledgling | 0.35 |
Using the default values here:
\[ F = 1.8 \times 3.5 \times 0.5 \times 0.35 = 1.1025 \]
This means each adult female is expected to contribute about 1.1 female fledglings per year in this simplified model.
Each factor represents one component of annual reproductive output, and together they combine multiplicatively to determine expected fecundity. For example, if a female attempts 1.8 nests per season and each nest has a 35% chance of success, then the expected number of successful nests is 0.63. Expected annual production of female fledglings then depends further on clutch size and the proportion of offspring that are female.
If a female attempts more nests, lays more eggs per attempt, and a greater fraction of nests succeed, her expected annual contribution of female offspring increases through the product of those components.
We will start by creating a custom function that we can use to construct a projection matrix:
build_matrix <-
function(
breeding_attempts,
clutch_size,
sex_ratio,
nest_success,
phi_j,
phi_a
) {
# Calculate fecundity:
fecundity <- breeding_attempts * clutch_size * sex_ratio * nest_success
# Create and return matrix:
c(
0,
phi_j,
fecundity,
phi_a
) %>%
matrix(
nrow = 2,
ncol = 2,
dimnames =
list(
c("Juvenile", "Adult"),
c("Juvenile", "Adult")
)
)
}… and then run our function based on our inputs:
projection_mat <-
build_matrix(
breeding_attempts = v_rate$breeding_attempts,
clutch_size = v_rate$clutch_size,
sex_ratio = v_rate$sex_ratio,
nest_success = v_rate$nest_success,
phi_j = v_rate$phi_j,
phi_a = v_rate$phi_a
)Let’s have a look:
| Juvenile | Adult | |
|---|---|---|
| Juvenile | 0.00 | 1.102 |
| Adult | 0.25 | 0.550 |
I know. The goal here is not to master linear algebra vocabulary (try not to fixate on it!), but to understand what these outputs tell us about population behavior.
To understand how the population changes over time, we perform an eigenanalysis — a mathematical procedure that extracts special numbers (eigenvalues) and vectors (eigenvectors) that describe the long-term behavior of the matrix.
Imagine using these vital rates to project the population forward many years. At first, the proportion of juveniles and adults may fluctuate depending on the starting population. But eventually the population settles into a stable mix of stages. This is called the stable stage distribution. Once that happens, the entire population simply grows or declines by the same factor every year. That factor is the population growth rate lambda \(\lambda\) and is defined mathematically as the dominant eigenvalue (i.e., the most biologically meaningful eigenvalue).
Lambda (\(\lambda\)) is the long-run annual population growth rate after the stage structure stabilizes. It tells you the multiplicative factor by which the population changes each year once it reaches a stable stage distribution.
| Value | Interpretation |
|---|---|
| \(\lambda > 1\) | Population is growing |
| \(\lambda = 1\) | Population is stable (replacement-level reproduction) |
| \(\lambda < 1\) | Population is declining |
For example, \(\lambda = 0.95\) means the population declines by 5% per year.
An eigenvector is a special direction in mathematical space that does not rotate when you apply the matrix — it only scales. In demographic terms:
Reproductive value is not the same as current reproduction. It reflects how much future population growth an individual in that stage is expected to generate on average.
Higher reproductive value means that individuals in that stage contribute more strongly to future population growth on average. Adults generally have higher reproductive value than juveniles because they have already survived the risky early-life period and are immediately available to breed.
We can break this into four steps:
Here we go!
# Get eigenvalues and eigenvectors of the projection matrix:
ev <- eigen(projection_mat)# Determine the id with the maximum eigenvalue:
lambda_idx <-
ev$values %>%
# `Re()` is used to ensure real values are returned (because `eigen()` can
# sometimes return complex numbers):
Re() %>%
# Get the index at the maximum eigenvalue:
which.max()
# Calculate lambda (population growth rate):
lambda <-
ev$values[lambda_idx] %>%
Re()Note: This sums to 1 (so it reads as proportions).
# Right eigenvector: stable stage distribution:
w <-
# Select the eigenvector associated with lambda:
ev$vectors[, lambda_idx] %>%
Re() %>%
# Normalize the vector to sum to 1:
{. / sum(.)}This uses a standard scaling convention (\(\mathbf{v} \cdot \mathbf{w} = 1\)) that ensures sensitivity and elasticity formulas work correctly. This is known as a dot product – a dot product is the sum of pairwise products across two vectors, written as \(\sum_i x_i y_i\).
# Left eigenvector: reproductive values:
ev_left <-
projection_mat %>%
# Transpose matrix (flip the rows and columns):
t() %>%
# Compute eigenvalues and eigenvectors:
eigen()
v <-
ev_left$vectors[, lambda_idx] %>%
Re() %>%
{. / as.numeric(. %*% w)}Consider the vector x = (2, 3) and y = (4, 5).
The dot product is:
\(\mathbf{x} \cdot \mathbf{y} = 2\times4 + 3\times5 = 23\)
| Quantity | Value |
|---|---|
| Growth rate | |
| $$\lambda$$ | 0.868 |
| $$log(\lambda) = r$$ | -0.142 |
| Stable stage distribution (w) | |
| Stable stage: Juvenile | 0.560 |
| Stable stage: Adult | 0.440 |
| Reproductive values (v) | |
| Reproductive value: Juvenile | 0.479 |
| Reproductive value: Adult | 1.662 |
Interpreting this output: A \(\lambda\) of 0.868 means the modeled population is declining. Under the stable stage distribution, roughly 56% of individuals are juveniles and 44% are adults. Adults have a higher reproductive value than juveniles, which makes biological sense: an adult has already cleared the high-mortality post-fledgling period and is positioned to breed immediately.
Lambda tells us whether the population grows or declines, while the eigenvectors describe the long-term stage structure and stage importance.
I know that formulas can be scary. On a first read, it is fine to focus on the biological meaning and the code comments, then come back to the notation after the big picture makes sense.
The formula below answers the question: if we modify a matrix element, or the vital rate that enters it, by a small absolute amount, how much would \(\lambda\) change?
This is known as Sensitivity – formally, the sensitivity of \(\lambda\) to matrix element \(a_{ij}\) is:
\[ s_{ij} = \frac{\partial \lambda}{\partial a_{ij}} = \frac{v_i \, w_j}{\mathbf{v} \cdot \mathbf{w}} \]
This is a product of reproductive values (how strongly individuals in each stage contribute to future population growth) and stable stage proportion (how common that stage is), normalized so the values are comparable. A high sensitivity means that even a small absolute nudge to that matrix element would substantially change \(\lambda\). Sensitivity is high when a stage is both common in the stable population and important for future growth.
In code, this formula is implemented by taking the outer product of reproductive values and stable stage proportions, then dividing by their dot product for normalization. In R, this is expressed as (v %o% w) / as.numeric(v %*% w) where:
v %o% w: The product of every possible combination of v and w (outer product of vectors). Note: An outer product of vectors is a matrix representing all pairwise products across two input vectors.as.numeric(v %*% w): The dot product of v and w, returned as a single numeric value for normalization.Consider the vector x = (2, 3) and y = (4, 5).
The outer product of vectors is:
\[ x_i \, y_j = \begin{pmatrix} 2\times4 & 2\times5 \\ 3\times4 & 3\times5 \end{pmatrix} = \begin{pmatrix} 8 & 10 \\ 12 & 15 \end{pmatrix} \]
Together, these two steps recreate the sensitivity formula: the outer product gives the numerator, and the dot product gives the denominator used for scaling.
# Sensitivity analysis:
(v %o% w) / as.numeric(v %*% w)Vital rates are measured on very different scales — a survival probability ranges from 0 to 1, while fecundity can easily exceed 2 or 3. Comparing raw sensitivities across vital rates that operate on different scales is misleading, like comparing miles per hour to kilometers per hour without converting. That is why we use elasticity instead.
The formula below answers the question: if a vital rate changes by 1%, how much would \(\lambda\) change in percentage terms?
This is known as Elasticity, a proportional (scale-free) version of sensitivity (Note: Elasticities sum to 1) that is formally described as:
\[ e_{ij} = \frac{a_{ij}}{\lambda} \cdot s_{ij} \]
Because elasticities are dimensionless (proportion in, proportion out), they are directly comparable across vital rates. A key property: the elasticities of all matrix elements sum to exactly 1.0, which means they divide the total proportional influence on \(\lambda\) among the matrix elements.
Sensitivity and elasticity are related, but they are not the same thing.
For example:
Suppose a vital rate changes from 0.50 to 0.51. That is an absolute change of 0.01, but a proportional change of 2%. Sensitivity uses the 0.01 framing and elasticity uses the 2% framing.
That is why elasticity is usually better for comparing very different kinds of vital rates.
Also note: A vital rate can have a high elasticity even if its current value is not large. Elasticity is about influence on lambda, not about which rate has the biggest numeric value.
Nest success (\(NS\)) enters the matrix through fecundity (\(F = B\times CS \times SR \times NS\)). Because all four components multiply together into a single matrix element, the elasticity of \(\lambda\) to each component individually equals the elasticity of \(\lambda\) to the whole fecundity element. This follows from the chain rule of calculus, but the important idea is simpler – because fecundity is built by multiplying several components together, a proportional change in any one component creates the same proportional change in fecundity.
The chain rule of calculus states that the proportional sensitivity of a product to any one of its multiplicative components is the same, regardless of that component’s scale.
# First, calculate the sensitivity matrix from the quantities that we already
# estimate above:
{(v %o% w) / as.numeric(v %*% w)} %>%
list(
# Then rescale the sensitivity matrix to generate the elasticity matrix:
(projection_mat / lambda) * .
) %>%
set_names("sensitivity_matrix", "elasticity_matrix") %>%
# Add column and row names:
map(
\(mat) {
dimnames(mat) <- dimnames(projection_mat)
mat
}
) %>%
list2env(.GlobalEnv)
# Extract values from the elasticity matrix:
elasticity_fecundity <- elasticity_matrix[1, 2]
elasticity_phi_j <- elasticity_matrix[2, 1]
elasticity_phi_a <- elasticity_matrix[2, 2]
# Put them together:
v_rate_elasticities <-
tibble(
vital_rate =
c(
"Nest success",
"Post-fledgling survival",
"Adult survival"
),
parameter =
c(
"nest_success",
"phi_j",
"phi_a"
),
value =
c(
v_rate$nest_success,
v_rate$phi_j,
v_rate$phi_a
),
elasticity =
c(
elasticity_fecundity,
elasticity_phi_j,
elasticity_phi_a
),
sensitivity =
c(
lambda * elasticity_fecundity / v_rate$nest_success,
lambda * elasticity_phi_j / v_rate$phi_j,
lambda * elasticity_phi_a / v_rate$phi_a
)
)| Rank | Vital rate | Parameter | Current value | Sensitivity | Elasticity | % of total |
|---|---|---|---|---|---|---|
| 1 | Adult survival | phi_a | 0.55 | 0.7320 | 0.464 | 46.4% |
| 2 | Nest success | nest_success | 0.35 | 0.6644 | 0.268 | 26.8% |
| 3 | Post-fledgling survival | phi_j | 0.25 | 0.9301 | 0.268 | 26.8% |
We want to know which biological rate matters most for population growth.
If we slightly improve a vital rate:
Elasticity is usually what we care about because it lets us compare very different kinds of rates (survival vs reproduction). For example:
If adult survival increases from 0.55 → 0.56 that is a +0.01 absolute change.
Sensitivity asks: “How much does \(\lambda\) change if survival increases by 0.01?”
Elasticity asks: “How much does \(\lambda\) change if survival increases by 1%?”
Because elasticity uses percentages, we can compare it fairly to fecundity or nest success.
Elasticity is usually more useful than sensitivity because it lets us compare very different vital rates on the same proportional scale.
Elasticity values describe the local slope of \(\lambda\) at the current vital rate estimates. But what does \(\lambda\) look like across the full plausible range of each vital rate? Perturbation analysis answers this by systematically varying one vital rate at a time while holding the others fixed at their baseline values, then recalculating \(\lambda\) at each step.
We’ll start by creating a function for calculating lambda from a population matrix:
get_lambda <-
function(mat) {
ev <- eigen(mat)
max(
Re(ev$values)
)
}Next, we’ll write a function that will get lambda across all potential values of a given vital rate:
perturb_grid <-
function(
v_rate_name,
seq_vals,
base_v_rate
) {
map_dfr(
seq_vals,
\(x) {
new_v_rate <- base_v_rate
new_v_rate[[v_rate_name]] <- x
mat <-
build_matrix(
breeding_attempts = new_v_rate$breeding_attempts,
clutch_size = new_v_rate$clutch_size,
sex_ratio = new_v_rate$sex_ratio,
nest_success = new_v_rate$nest_success,
phi_j = new_v_rate$phi_j,
phi_a = new_v_rate$phi_a
)
tibble(
vital_rate = v_rate_name,
x = x,
lambda = get_lambda(mat)
)
}
)
}… and then calculate using map_df() to iterate across potential values:
perturbation_data <-
c(
"nest_success",
"phi_j",
"phi_a"
) %>%
map_df(
~ perturb_grid(
.x,
seq(
0.05,
0.95,
by = 0.01
),
v_rate
)
) %>%
mutate(
vital_rate_label =
case_when(
vital_rate == "nest_success" ~ "Nest success",
vital_rate == "phi_j" ~ "Post-fledgling survival",
vital_rate == "phi_a" ~ "Adult survival"
) %>%
fct_relevel(
"Adult survival",
"Nest success",
"Post-fledgling survival"
)
)… and get the points at our current vital rate values:
baseline_pts <-
tibble(
vital_rate_label =
c(
"Nest success",
"Post-fledgling survival",
"Adult survival"
) %>%
fct_relevel(
"Adult survival",
"Nest success",
"Post-fledgling survival"
),
x =
c(
v_rate$nest_success,
v_rate$phi_j,
v_rate$phi_a
),
lambda = lambda
)v_rate_elasticities %>%
mutate(
vital_rate =
fct_reorder(
vital_rate,
elasticity,
.desc = TRUE
)
) %>%
ggplot() +
aes(
x = vital_rate,
y = elasticity,
fill = vital_rate
) +
geom_col(
width = 0.55,
show.legend = FALSE
) +
geom_text(
aes(
label = glue::glue("e = {round(elasticity, 3)}")
),
vjust = -0.4,
size = 3.8
) +
scale_fill_manual(
values =
c(
"Adult survival" = "#2166AC",
"Nest success" = "#D6604D",
"Post-fledgling survival" = "#4DAC26"
)
) +
scale_y_continuous(
limits = c(0, 0.80),
expand = c(0, 0),
labels =
label_number(accuracy = 0.01)
) +
labs(
title = "Vital-rate elasticities",
subtitle =
glue("\u03bb = {round(lambda, 3)} | Post-breeding census, 2-stage model"),
x = NULL,
y = "Elasticity of \u03bb"
) +
theme_bw(base_size = 13) +
theme(
plot.title =
element_text(face = "bold"),
plot.subtitle =
element_text(color = "grey40")
)
How to read this chart: Each bar represents one focal vital rate. The height of the bar is the elasticity — the proportional change in \(\lambda\) you would expect from a 1% proportional change in that vital rate. Because elasticities sum to 1, the bars can be interpreted like slices of a pie: a bar at 0.55 means that vital rate accounts for 55% of the total proportional influence represented in this model.
The vital rate with the tallest bar is the one that, proportionally speaking, exerts the strongest proportional influence on population growth within the modeled life cycle.
If the elasticity of adult survival is 0.46, this means that a 1% increase in adult survival would increase λ by approximately 0.46%.
Likewise, if the elasticity of nest success is 0.27, a 1% increase in nest success would increase \(\lambda\) by about 0.27%.
So a proportional change in adult survival would have about 1.73 times the effect on population growth as the same proportional change in nest success.
As you interpret the curves below, ask yourself three questions:
perturbation_data %>%
ggplot() +
aes(
x = x,
y = lambda,
color = vital_rate_label
) +
geom_line(linewidth = 1.1) +
geom_hline(
yintercept = 1,
linetype = "dashed",
color = "grey50"
) +
geom_point(
data = baseline_pts,
aes(
x = x,
y = lambda
),
shape = 21,
size = 3,
fill = "white",
stroke = 1.5
) +
facet_wrap(~ vital_rate_label, scales = "free_x") +
scale_color_manual(
values =
c(
"Adult survival" = "#2166AC",
"Nest success" = "#D6604D",
"Post-fledgling survival" = "#4DAC26"
),
guide = "none"
) +
scale_y_continuous(
labels = label_number(accuracy = 0.01)
) +
labs(
title = "Lambda response to vital rate perturbation",
subtitle = "Open circle = baseline estimate | Dashed line = \u03bb = 1.0",
x = "Vital rate value",
y = "Lambda (\u03bb)"
) +
theme_bw(base_size = 13) +
theme(
plot.title =
element_text(face = "bold"),
plot.subtitle =
element_text(color = "grey40"),
strip.text =
element_text(face = "bold")
)
How to read these curves: Each panel shows a single vital rate on the x-axis and the resulting \(\lambda\) on the y-axis. The open circle marks the current (baseline) estimate. The dashed horizontal line at \(\lambda = 1\) marks the threshold between growth and decline.
Why this figure is useful: Elasticity tells us what happens for very small changes near the current estimate. The perturbation curves show what happens across the entire plausible range of the vital rate. This helps answer questions like: How much would adult survival need to increase before the population shifts from decline to growth?
Several things to look for:
| Rank | Vital rate | Current value | Elasticity | Sensitivity | % of total |
|---|---|---|---|---|---|
| 1 | Adult survival | 0.55 | 0.464 | 0.7320 | 46.4% |
| 2 | Nest success | 0.35 | 0.268 | 0.6644 | 26.8% |
| 3 | Post-fledgling survival | 0.25 | 0.268 | 0.9301 | 26.8% |
How to read this table:
This table identifies which demographic process exerts the greatest proportional influence on \(\lambda\) under the assumptions of the model. The top-ranked vital rate has the greatest leverage on population growth.
In this model, we can see that Adult survival has the greatest proportional influence on \(\lambda\).
A useful extension of our perturbation analysis is to solve directly for the value of a vital rate that would place the population exactly at replacement.
For example, we can use it to address the question ‘Given fixed adult and post-fledgling survival, what value of nest success would bring the population to exactly \(\lambda = 1\)?’ This can be solved analytically from the characteristic equation of the 2×2 matrix. Setting \(\lambda = 1\) and rearranging for \(NS\):
\[ NS^* = \frac{1 - \phi_a}{\phi_j \cdot B \cdot CS \cdot SR} \]
ns_required <-
(1 - v_rate$phi_a) /
(
v_rate$phi_j *
v_rate$breeding_attempts *
v_rate$clutch_size *
v_rate$sex_ratio
)
ns_change <- ns_required - v_rate$nest_successUnder the current vital rates, nest success would need to reach 0.571 — a change of 0.221 from its baseline value of 0.35 — for the population to achieve replacement-level growth. This represents an increase relative to the current estimate.
For a 2×2 matrix, setting \(\lambda = 1\) in the characteristic equation reduces to a linear equation in \(NS\), so the solution is exact. No iteration required. This approach works whenever \(NS\) enters the matrix through a single multiplicative fecundity term, as it does here.
In this tutorial, we built a simple two-stage population model, used eigenanalysis to describe long-term population behavior, and applied sensitivity and elasticity analysis to identify which demographic processes matter most for population growth. The key idea is that matrix models let us connect biological rates such as survival and nest success to overall population dynamics in a rigorous and interpretable way.
| Term | Definition |
|---|---|
| Adult survival (\(\phi_a\)) | Probability a breeding adult survives from one breeding season to the next |
| Chain rule | A calculus rule that describes how to differentiate a composite function; here it explains why each multiplicative component of fecundity shares the same elasticity |
| Dot product | Sum of the pairwise products across two vectors, written as \(\sum_i x_i y_i\) and expressed in R as as.numeric(x %*% y) |
| Elasticity | The proportional change in \(\lambda\) resulting from a proportional change in a vital rate; scale-free and sums to 1 across all matrix elements |
| Eigenvalue | A scalar (number) associated with a matrix that describes how much a special vector is stretched or compressed by the matrix; in demography, the dominant eigenvalue is \(\lambda\) |
| Eigenvector | A vector whose direction is unchanged by the matrix; in demography, the right eigenvector gives the stable stage distribution and the left eigenvector gives reproductive values |
| Fecundity (\(F\)) | Per-capita contribution of adults to the juvenile stage class; the product of breeding attempts, clutch size, sex ratio, and nest success |
| Lambda (\(\lambda\)) | Asymptotic annual population growth rate; the factor by which population size changes per year once stage proportions stabilize |
| log(\(\lambda\)) = \(r\) | The continuous-time intrinsic rate of increase; equivalent to lambda on a log scale |
| Matrix model | A mathematical model that projects population dynamics by multiplying a matrix of vital rates by a vector of stage-specific abundances |
| Nest success | Probability a nest produces ≥ 1 fledgling |
| Outer product of vectors | The product of every possible pair of elements from two input vectors, expressed in R as x %o% y |
| Perturbation analysis | A systematic examination of how \(\lambda\) responds to changes in a vital rate across its full plausible range |
| Post-breeding census | A census taken after the breeding season, so newly fledged young are included in the count |
| Post-fledgling survival (\(\phi_j\)) | Probability a fledgling survives from leaving the nest through its first winter to become a breeding adult |
| Reproductive value | The relative expected contribution of an individual in a given stage class to future population growth |
| Stable stage distribution | The equilibrium proportion of individuals in each stage class; the distribution the population converges to regardless of its starting composition |
| Sensitivity | The absolute change in \(\lambda\) resulting from a small absolute change in a vital rate or matrix element |