🧪 Part 2: Estimating Quasi-Extinction Risk in a Variable Environment

In this exercise, we’ll use matrix population models (MPMs) again to simulate how a population behaves under fluctuating environmental conditions. This time, we’ll focus on quasi-extinction risk — the chance that a population falls below a critical threshold, even if it hasn’t gone extinct.

👉 Quasi-extinction doesn’t mean the population has hit zero — just that it has dropped so low that recovery might be unlikely.

📦 Setup: Re-use the Code from Before

We’re going to use the same setup as in Part 1. If you’re running this code fresh, you’ll need to load everything again.

# Set seed for reproducibility
set.seed(42)

# Define projection matrices for Good and Bad years
A_good <- matrix(c(
  0.0, 1.2, 1.5,
  0.5, 0.3, 0.0,
  0.0, 0.6, 0.7
), nrow = 3, byrow = TRUE)

A_bad <- matrix(c(
  0.0, 0.8, 0.75,
  0.3, 0.2, 0.0,
  0.0, 0.4, 0.5
), nrow = 3, byrow = TRUE)

# Initial population: 50 juveniles, 30 subadults, 20 adults
n0 <- c(50, 30, 20)

# Simulation parameters
timesteps <- 100     # Years to simulate
n_sim <- 1000        # Number of simulations
p_good <- 0.5        # Probability that any year is a good year

# Set up storage for results
lambda_s <- numeric(n_sim)
trajectory_data <- matrix(NA, nrow = timesteps + 1, ncol = n_sim)

🔁 Run the Simulations Again

We’ll simulate 1000 possible futures for the population, switching randomly between good and bad years.

for (i in 1:n_sim) {
  pop <- n0
  total_N <- numeric(timesteps + 1)
  total_N[1] <- sum(pop)
  
  for (t in 1:timesteps) {
    A_t <- if (runif(1) < p_good) A_good else A_bad
    pop <- A_t %*% pop
    total_N[t + 1] <- sum(pop)
  }
  
  time <- 0:timesteps
  lambda_s[i] <- coef(lm(log(total_N) ~ time))[2]
  trajectory_data[, i] <- total_N
}

📉 Visualize the Results

Let’s plot all the simulation runs on a log scale to show variation in population size over time.

# Log transform the data
trajectory_data <- log(trajectory_data)

# Plot all runs
matplot(0:timesteps, trajectory_data, type = "l", lty = 1,
        col = adjustcolor("darkblue", alpha.f = 0.3),
        xlab = "Time (years)", ylab = "Total Population Size (log)",
        main = "Variation in Stochastic Trajectories (1000 Runs)")

# Add mean trajectory
lines(0:timesteps, rowMeans(trajectory_data, na.rm = TRUE), 
      col = "red", lwd = 2)

legend("topright", legend = c("Individual runs", "Mean"), 
       col = c("darkblue", "red"), lty = c(1, 1), lwd = c(1, 2), bty = "n")

❗ Quasi-Extinction Risk Analysis

Here’s the new part! We’re going to look at how many simulations drop below a critical population size — our quasi-extinction threshold.

Step 1: Choose a time frame

We’ll check for quasi-extinction within the first 50 years.

timeLimit <- 50
trajectory_data_cut <- trajectory_data[1:timeLimit, ]

Step 2: Convert population sizes back from log-scale

exp_data <- exp(trajectory_data_cut)

Step 3: Set a threshold

Let’s say a population size below 10 individuals counts as quasi-extinct.

threshold <- 10

Step 4: Check which simulations dipped below the threshold

We’ll look across all 50 years for each simulation to see if the population ever dropped below 10.

quasiExtinct <- apply(exp_data, 2, function(x) any(x < threshold))

# Count how many simulations show quasi-extinction
table(quasiExtinct)

## quasiExtinct
## FALSE  TRUE 
##   990    10

# As a proportion
table(quasiExtinct)[2] / sum(table(quasiExtinct))

## TRUE 
## 0.01

🔍 Interpretation

The proportion of simulations that dip below the threshold gives us a risk estimate. If many of them do, it may suggest a significant risk of collapse, even if the mean population looks stable.

🎛️ Things to Explore

Here are some what-if questions to try:

1. What happens if bad years become more common?

Try setting p_good <- 0.3 instead of 0.5.
This simulates climate change or deteriorating conditions.

2. What if you change the extinction threshold?

Try setting threshold <- 40 or even 60.
This simulates different definitions of what counts as a dangerously small population.

3. What if you extend or shorten the time window?

Try changing timeLimit <- 30 or timeLimit <- 100.
Populations might survive short-term, but collapse later.

4. Time to extinction Instead of just asking if extinction happens, record when it happens.

Use the following code, above:

time_to_extinction <- apply(exp_data, 2, function(x) {
  w <- which(x < threshold)
  if (length(w) == 0) return(NA)
  return(min(w))
})
hist(time_to_extinction, main = "Time to Quasi-Extinction", xlab = "Years")

5. Rescue mechanism (advanced!)

Add a “rescue” mechanism where, if the population falls below a certain level, it receives a boost (e.g., through migration or conservation intervention). This code should help.

total_N[t + 1] <- sum(pop) #find this line of code in your script
if(total_N[t + 1] < 100){pop <- pop + c(0,0,10)} #add this line

6. Use real data (very advanced!)

Obtain real world data from the COMADRE or COMPADRE databases and use those to run extinction risk estimates.
Use package Rcompadre and modify the code above to take an arbitrary number of matrices.

🧠 Key Takeaways

Quasi-extinction lets us assess risk before actual extinction occurs.
Environmental variability can lead to large differences between simulations.
Extinction risk is not simple — it must be carefully defined based on specific criteria, including both time (how far into the future we look) and numbers (what population size counts as “extinct”).
Even if the average trend looks okay, some populations may still crash — this is why stochastic simulations are so powerful!

Stochastic Matrix Models: Part 2

Owen Jones