Aggregating Asian Elephant Matrix Population Models for Comparative Demography with mpmaggregate

Richard A. Hinrichsen

ORCID: https://orcid.org/0000-0003-0761-3005

Overview

Age-structured matrix population models (MPMs) can differ in age-class width and projection interval, even when describing the same species. This creates a practical problem in comparative demography: how can models that represent life cycles at different temporal resolution be compared?

In this vignette, we demonstrate how to use the package mpmaggregate to aggregate a fine-grained, annual (1-year) Asian elephant Leslie MPM into a coarser, 5-year Leslie MPM, and then compare the aggregated model to an independently parameterized coarser Leslie MPM that represents a different population. Because age classes in a Leslie model are defined by the model’s time step, aggregating from 1-year to 5-year intervals necessarily redefines the age-class structure (i.e., age classes become 5-year bins) while producing transitions appropriate to the longer projection interval.

We illustrate this approach using two Asian elephant populations from COMADRE database (Salguero‐Gómez et al. 2016) accessed via Rcompadre (Jones et al. 2022):

Both MPMs are Leslie matrix models parameterized using a post-breeding census formulation. Our focus is the aggregation of the Leslie MPM for the Priiyar Reserve population with an annual projection interval to a Leslie MPM with a 5-year projection interval. To do this, we use leslie_aggregate(), which performs Leslie-to-Leslie MPM aggregation under four alternative assumptions defined by the choice of demographic framework ("lambda" or "R0") and criterion ("standard" or "elasticity"):

  1. framework = "lambda", criterion = "standard"
  2. framework = "lambda", criterion = "elasticity"
  3. framework = "R0", criterion = "standard"
  4. framework = "R0", criterion = "elasticity"

criterion = "elasticity" implies elasticity-consitent aggregation (Hinrichsen et al., 2026).

Throughout, we seek to uncover the substantive demographic differences between the two Asian elephant populations, using \(\lambda\), \(R_0\), generation times, and elasticities. We also examine at how results differ depending on framework ("lambda" or "R0") and criterion ("standard" or "elasticity").

Loading packages

We begin by loading the mpmaggregate package, which provides the functions used throughout this vignette for aggregating matrix population models under the \(\lambda\) and \(R_0\) demographic frameworks. The package knitr is used for report generation, and kableExtra is used for table formatting.

library(mpmaggregate)
library(kableExtra)
library(rphylopic)
#> You are using rphylopic v.1.6.0. Please remember to credit PhyloPic contributors (hint: `get_attribution()`) and cite rphylopic in your work (hint: `citation("rphylopic")`).

Data acquisition (COMADRE via Rcompadre)

The Asian Elephant data is stored in the data directory. Retrieve data using data() command.


# Matrices for Asian elephant populations included in this package were
# originally retrieved from COMADRE using the Rcompadre package.

# Not run: requires Rcompadre and internet access

# library(Rcompadre)
# comadre <- cdb_fetch("comadre",version = "4.23.3.1")
# mpm_elephant1 <- comadre[comadre$MatrixID == 249273, ]
# matA_elephant1 <- matA(mpm_elephant1)[[1]]             # Projection interval is 5 years
# matA_elephant1 <- matA_elephant1[1:12,1:12]            # Confine to females 
# mpm_elephant2 <- comadre[comadre$MatrixID == 249274, ]
# matA_elephant2 <- matA(mpm_elephant2)[[1]]             # Projection interval is 1 year
# matA_elephant2[60, 60] <- 0                            # Remove statis to form a true Leslie Matrix
# usethis::use_data(matA_elephant1, matA_elephant2,      # Write data (development)
# overwrite = TRUE)    

# Load data locally
data(matA_elephant1, matA_elephant2, package="mpmaggregate")
matA_coarse <- matA_elephant1 # 12x12 5-year projection interval
matA_fine <- matA_elephant2   # 60x60 1-year projection interval

Define the aggregation target

We aggregate the Periyar Reserve Asian elephant MPM from an annual time step to a 5-year time step, allowing us to directly compare its vital rates to those of the Nagarhole National Park Asian elephant MPM.

Rigorous methods for Leslie-to-Leslie MPM aggregation are implemented using the function leslie_aggregate(), which provides four aggregation options. If the goal is to preserve the asymptotic population growth rate, one uses framework = "lambda". Alternatively, if the goal is to preserve the net reproductive rate, one uses framework = "R0".

Within either framework, aggregation is performed using a standard or an elasticity-consistent criterion. The standard aggregator ensures consistency of the stable age distribution implied by the model, while the elasticity-consistent aggregator additionally ensures consistency of reproductive values. In the R0 framework, the stable age distribution is the cohort stable age distribution, and reproductive values are the cohort reproductive values.

We now run leslie_aggregate() under the four combinations.

Run aggregation under four option sets


res <- leslie_aggregate(
  matA_fine,
  ngroups = 12,
  framework = "lambda",
  criterion = "standard"
)
matB_lambda_stand <- res$matAk_agg
matB_lambda_stand_eff <- res$effectiveness

res <- leslie_aggregate(
  matA_fine,
  ngroups = 12,
  framework = "lambda",
  criterion = "elasticity"
)
matB_lambda_elast <- res$matAk_agg
matB_lambda_elast_eff <- res$effectiveness

res <- leslie_aggregate(
  matA_fine,
  ngroups = 12,
  framework = "R0",
  criterion = "standard"
)
matB_R0_stand <- res$matAk_agg
matB_R0_stand_eff<- res$effectiveness

res <- leslie_aggregate(
  matA_fine,
  ngroups = 12,
  framework = "R0",
  criterion = "elasticity"
)
matB_R0_elast <- res$matAk_agg
matB_R0_elast_eff <- res$effectiveness

Comparison to an existing coarse (5-year) matrix

The COMADRE matrix for Asian elephants in the Nagarhole National Park already has ProjectionInterval == 5, so we compare the aggregated 5-year matrices for Asian elephants in the Periyar Reserve against it.

Demographic metrics for Leslie matrices

For Leslie matrices, fertility rates occupy the first row, while survival probabilities occupy the subdiagonal. This structure allows direct decomposition of the projection matrix \(\mathbf{A}\) into survival probability \(\mathbf{U}\)and fertility rate \(\mathbf{F}\)components.


# Decompose a Leslie matrix into its components U and F
leslie_UF <- function(A) {
  F <- matrix(0, nrow(A), ncol(A))
  F[1, ] <- A[1, ]
  U <- A
  U[1, ] <- 0
  list(U = U, F = F)
}

# Return R0, the net reproductive rate
get_R0_leslie <- function(A) {
  uf <- leslie_UF(A)
  U <- uf$U
  F <- uf$F
  
  I <- diag(nrow(U))
  N <- solve(I - U)
  K <- F %*% N
  spectral_radius(K)
}

# Return generation time in the lambda framework
get_generation_time_lambda <- function(A) {
  uf <- leslie_UF(A)
  U <- uf$U
  F <- uf$F
  I <- diag(nrow(U))
  N <- solve(I - U)
  gt = generation_time(F, U, framework = "lambda")$generation_time
  return(gt)
}

# Return cohort generation time
get_generation_time_R0 <- function(A) {
  uf <- leslie_UF(A)
  U <- uf$U
  F <- uf$F
  I <- diag(nrow(U))
  N <- solve(I - U)
  gt = generation_time(F, U, framework = "R0")$generation_time
  return(gt)
}

Comparison of λ and \(R_0\) across five models. First develop a list of models with associated projection matrices

models <- list(
  "Orig. coarse"             = matA_coarse,
  "Agg λ + standard"         = matB_lambda_stand,
  "Agg λ + elasticity"       = matB_lambda_elast,
  "Agg R0 + standard"        = matB_R0_stand,
  "Agg R0 + elasticity"      = matB_R0_elast
)

Next table the results for λ and \(R_0\)


#Form the data frame that contains the table information
results <- data.frame(
  Model  = names(models),
  lambda = sapply(models, spectral_radius),
  R0     = sapply(models, get_R0_leslie),
  Ta     = sapply(models, get_generation_time_lambda)*5,  #measure in years
  Tc     = sapply(models, get_generation_time_R0)*5,      #measure in years
  Eff    = c(NA,matB_lambda_stand_eff,matB_lambda_elast_eff,
             matB_R0_stand_eff,matB_R0_elast_eff),
  row.names = NULL,
  stringsAsFactors = FALSE
)

n_rows <- nrow(results)

#Construct a table using function kable
knitr::kable(
  results,
  col.names = c(
    "Model",
    "<em>&lambda;</em>",
    "<em>R</em><sub>0</sub>",
    "<em>T</em><sub>a</sub> (years)",
    "<em>T</em><sub>c</sub> (years)",
    "<em>&rho;</em><sup>2</sup>"
  ),
  digits = 3,
  caption = paste0(
  "<strong>Table 1.</strong> Demographic comparison between two Asian elephant populations, ",
  "represented by an original coarse-stage Leslie model from Nagarhole National Park ",
  "(MatrixID 249273) and aggregated versions of a Leslie model from Periyar Reserve ",
  "(MatrixID 249274). Metrics include population growth rate (<em>&lambda;</em>), ",
  "net reproductive rate (<em>R</em><sub>0</sub>), generation time (years; ",
  "<em>T</em><sub>a</sub>), cohort generation time (years; ",
  "<em>T</em><sub>c</sub>), and aggregation effectiveness (<em>&rho;</em><sup>2</sup>). ",
  "&ldquo;Standard&rdquo; denotes use of a standard aggregator and ",
  "&ldquo;elasticity&rdquo; denotes elasticity-consistent aggregation."
  ),
  format = "html",
  escape = FALSE
) |>
  kable_styling(full_width = TRUE) |>
  row_spec(0, extra_css = "border-bottom: 2px solid black;") |>
  row_spec(nrow(results), extra_css = "border-bottom: 2px solid black;") |>
  footnote(
    general_title = "",
    general = paste0(
      "<span style='font-size: 90%;'>",
      "<strong>Note.</strong> The projection interval associated with the population ",
      "growth rate (<em>&lambda;</em>) is 5 years.</span>"
    ),
    footnote_as_chunk = TRUE,
    escape = FALSE
  )
Table 1. Demographic comparison between two Asian elephant populations, represented by an original coarse-stage Leslie model from Nagarhole National Park (MatrixID 249273) and aggregated versions of a Leslie model from Periyar Reserve (MatrixID 249274). Metrics include population growth rate (λ), net reproductive rate (R0), generation time (years; Ta), cohort generation time (years; Tc), and aggregation effectiveness (ρ2). “Standard” denotes use of a standard aggregator and “elasticity” denotes elasticity-consistent aggregation.
Model λ R0 Ta (years) Tc (years) ρ2
Orig. coarse 1.114 2.157 34.055 37.279 NA
Agg λ + standard 1.291 4.784 27.596 34.123 0.990
Agg λ + elasticity 1.291 4.836 27.813 34.332 0.996
Agg R0 + standard 1.293 4.811 27.447 34.018 0.992
Agg R0 + elasticity 1.292 4.811 27.581 34.128 0.993
Note. The projection interval associated with the population growth rate (λ) is 5 years.

How to read Table 1. Across aggregation approaches, estimates of \(\lambda\) and \(R_0\) for the Periyar Reserve population are similar and indicate faster population growth and shorter generation time relative to the Nagarhole National Park population. Each row corresponds to a Leslie matrix population model: the original coarse model from Nagarhole National Park (COMADRE) and four aggregated versions derived from the annual Periyar Reserve model. Differences among the aggregated models reflect the aggregation objective (preserving \(\lambda\) versus preserving \(R_0\)) and the fitting criterion (standard versus elasticity-consistent). Values of aggregation effectiveness, \(\rho^2\), indicate good agreement with the original model for the Periyar Reserve population.

Interpretation notes

Plot figures

Figure 1. Vital rates across five Leslie models. Panel (a) shows fertility rates (first-row entries of the Leslie matrix, \(F_i\)). Panel (b) shows survival probabilities (subdiagonal entries, \(S_i = A_{i+1,i}\)). Bars are grouped by vital rate; within each group, colors correspond to the five different models (original coarse model for the Nagarhole National Park population
and four 5-year aggregations of the model for the Periyar Reserve population). Only vital rates that are nonzero in at least one model are shown. The models in the Figure Legend are as described in the caption of Table 1.


#Asian elephant
#"Elephas maximus" by  T. Michael Keesey (2011; CC BY 3.0).

# Not run: requires internet access

# uuid <- "7c9ab182-175d-4f02-96d0-09c1e5212bff"
# img <- get_phylopic(uuid = uuid, format = "raster", height = 691)
# Write image 
# png::writePNG(unclass(img), "inst/extdata/elephant.png")

#Load the image file that is stored locally 
img_file <- system.file("extdata", "elephant.png", package = "mpmaggregate")

if (!nzchar(img_file)) {
  img_file <- "../inst/extdata/elephant.png"
}

img <- png::readPNG(img_file)
# ---- Helpers to extract Leslie vital rates from A ----
extract_fertility <- function(A) {
  as.numeric(A[1, ])
}

extract_survival <- function(A) {
  n <- nrow(A)
  s_sub  <- A[cbind(2:n, 1:(n - 1))]  # A[i+1, i]
  c(s_sub)
}

model_names <- names(models)

# Choose distinct colors for the 5 models (base R friendly)
cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e")

# Build rate-by-model matrices so:
# rows = models (5 bars), cols = vital rates (groupings)
F_by_model <- t(sapply(models, extract_fertility))  # 5 x n
S_by_model <- t(sapply(models, extract_survival))   # 5 x n-1

# Name columns (vital rates)
n <- nrow(models[[1]])
colnames(F_by_model) <- paste0("F", 1:n)
colnames(S_by_model) <- c(paste0("S", 1:(n - 1)))

# Keep only vital rates that are nonzero in at least one model
keep_F <- colSums(abs(F_by_model) > 0) > 0
keep_S <- colSums(abs(S_by_model) > 0) > 0

F_plot <- F_by_model[, keep_F, drop = FALSE]
S_plot <- S_by_model[, keep_S, drop = FALSE]

# Leave extra space on the right for a legend that won't overlap bars
op <- par(
  mfrow = c(2, 1),
  mar = c(7, 4, 3, 7),
  oma = c(0, 0, 0, 0)
)

op <- par(
  mfrow = c(2, 1),
  mar = c(5, 4, 1.5, 7),
  oma = c(0, 0, 0, 0)
)

ylim = c(0, max(F_plot))
mydiff = ylim[2] - ylim[1]
ylim[2] = ylim[2] + .3 * mydiff

# Panel (a): Fertility rates
barplot(
  F_plot,
  ylim = ylim,
  beside = TRUE,
  col = cols,
  border = "white",
  las = 2,
  ylab = "Fertility rate",
  main = "Fertility rates"
)
mtext("(a)",
      side = 3,
      adj = 0,
      line = 0.2)
lims <-  par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x =  NULL,
  y = mydiff * .9,
  height =  .3 * mydiff * dim(img)[1] / 1024
)

# One legend for both panels, placed in the right margin of the top panel
par(xpd = NA)
legend(
  "topright",
  inset = c(-0.18, 0.00),
  legend = model_names,
  fill = cols,
  bty = "n",
  cex = 0.9
)
par(xpd = FALSE)
ylim <-  c(0, max(S_plot))
mydiff <- ylim[2] - ylim[1]
ylim[2] <- ylim[2] + .3 * mydiff

# Panel (b): Survival probs
barplot(
  S_plot,
  ylim = ylim,
  beside = TRUE,
  col = cols,
  border = "white",
  las = 2,
  ylab = "Survival probability",
  main = "Survival probabilities"
)
mtext("(b)",
      side = 3,
      adj = 0,
      line = 0.2)
lims <-  par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x =  NULL,
  y = mydiff * .9,
  height =  .3 * mydiff * dim(img)[1] / 1024
)

# Legend for bottom panel (same placement logic)
par(xpd = NA)
legend(
  "topright",
  inset = c(-0.18, 0.00),
  legend = model_names,
  fill = cols,
  bty = "n",
  cex = 0.9
)

par(xpd = FALSE)

par(op)

How to read Figure 1. Each x-axis label is a vital rate. Fertility rate, \(F_i\) is the contribution of age class \(i\) to newborns over the projection interval. Survival probability, \(S_i\) is the probability of surviving and advancing from age class \(i\) to \(i+1\) over the projection interval. Because the four aggregated matrices are all 5-year models, differences among them reflect the aggregation objective (\(\lambda\) vs \(R_0\)) and the fitting criterion (standard vs elasticity-consistent), rather than differences in projection interval. The figure shows that the Nagarhole National Park population has lower population growth rate than the Periyar Reserve population due to lower fertility rates.The silhouette is from PhyloPic (https://www.phylopic.org; Keesey, 2023) and were added using the rphylopic R package (Gearty & Jones, 2023). The silhouette was contributed as follows: “Elephas maximus” by T. Michael Keesey (2011; CC BY 3.0).

Figure 2. Elasticities of \(\lambda\) with respect to matrix entries. Elasticities quantify the proportional effect on \(\lambda\) of proportional changes in each vital rate. Panel (a) shows elasticities with respect to fertility rates (first row), and panel (b) shows elasticities with respect to survival probabilities (subdiagonal). Bars are grouped by vital rate, with colors indicating the five different models (original coarse model for the Nagarhole National Park population
and four 5-year aggregations of the model for the Periyar Reserve population). Elasticities sum to 1 within each matrix, so the panels show how the influence on \(\lambda\) is distributed across life-history components. The models in the figure legend are as described in the caption of Table 1.

#Asian elephant
#"Elephas maximus" by  T. Michael Keesey (2011; CC BY 3.0).

# Not run: requires internet access

# uuid <- "7c9ab182-175d-4f02-96d0-09c1e5212bff"
# img <- get_phylopic(uuid = uuid, format = "raster", height = 691)
# Write image 
# png::writePNG(unclass(img), "inst/extdata/elephant.png")

#Load a local copy of the image file
img_file <- system.file("extdata", "elephant.png", package = "mpmaggregate")

if (!nzchar(img_file)) {
  img_file <- "../inst/extdata/elephant.png"
}

img <- png::readPNG(img_file)

# Compute elasticity matrices for each model (same dimensions as A)
E_models <- lapply(models, function(A) {
  E <- mpm_elasticity(A, framework = "lambda")$elasticity
  as.matrix(E)
})

# ---- Helpers to extract Leslie elasticities from an elasticity matrix ----
extract_fertility_elast <- function(E) {
  as.numeric(E[1, ])
}

extract_survival_elast <- function(E) {
  n <- nrow(E)
  e_sub  <- E[cbind(2:n, 1:(n - 1))]  # elasticity of A[i+1, i]
  as.numeric(e_sub)
}

model_names <- names(models)

# Same colors as the vital-rate plot for consistency
cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e")

# Build elasticity-by-model matrices so:
# rows = models (5 bars), cols = vital rates (groupings)
F_el_by_model <- t(sapply(E_models, extract_fertility_elast))  # 5 x n
S_el_by_model <- t(sapply(E_models, extract_survival_elast))   # 5 x n-1

# Name columns (vital rates)
n <- nrow(models[[1]])
colnames(F_el_by_model) <- paste0("F", 1:n)
colnames(S_el_by_model) <- c(paste0("S", 1:(n - 1)))

# Keep only elasticities that are nonzero in at least one model
keep_F <- colSums(abs(F_el_by_model) > 0) > 0
keep_S <- colSums(abs(S_el_by_model) > 0) > 0

F_el_plot <- F_el_by_model[, keep_F, drop = FALSE]
S_el_plot <- S_el_by_model[, keep_S, drop = FALSE]

# Leave extra space on the right for legends that won't overlap bars
op <- par(
  mfrow = c(2, 1),
  mar = c(7, 4, 3, 7),
  oma = c(0, 0, 0, 0)
)

op <- par(
  mfrow = c(2, 1),
  mar = c(5, 4, 1.5, 7),
  oma = c(0, 0, 0, 0)
)

# Panel (a): Elasticities with respect to fertility rates
barplot(
  F_el_plot,
  beside = TRUE,
  ylim = c(0, 0.20),
  col = cols,
  border = "white",
  las = 2,
  ylab = "Elasticity of λ",
  main = "Elasticities with respect to fertility rates"
)
mtext("(a)",
      side = 3,
      adj = 0,
      line = 0.2)
lims <-  par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x =  NULL,
  y = mydiff * .75,
  height =  .3 * mydiff * dim(img)[1] / 1024
)
par(xpd = NA)
legend(
  "topright",
  inset = c(-0.15, -0.05),
  legend = model_names,
  fill = cols,
  bty = "n",
  cex = 0.9
)
par(xpd = FALSE)

# Panel (b): Elasticities with respect to survival probabilities
barplot(
  S_el_plot,
  beside = TRUE,
  ylim = c(0, .20),
  col = cols,
  border = "white",
  las = 2,
  ylab = "Elasticity of λ",
  main = "Elasticities with respect to survival probabilities"
)
mtext("(b)",
      side = 3,
      adj = 0,
      line = 0.2)
lims <-  par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x =  NULL,
  y = mydiff * .75,
  height =  .3 * mydiff * dim(img)[1] / 1024
)

par(xpd = NA)
legend(
  "topright",
  inset = c(-0.15, -0.05),
  legend = model_names,
  fill = cols,
  bty = "n",
  cex = 0.9
)

par(xpd = FALSE)

par(op)

How to read Figure 2. A larger elasticity (of lambda with respect to matrix entries) for a given vital rate means that small proportional changes to that vital rate would have a comparatively larger effect on long-term population growth. For Asian elephants, elasticities are concentrated in survival probabilities rather than fertility rates, reflecting the strong influence of survivorship on \(\lambda\). Aggregation methods yield similar elasticities. Fertility rates and late-life survival probabilities have higher elasticities in the Nagarhole National Park population than in the Periyar Reserve population.The silhouette is from PhyloPic (https://www.phylopic.org; Keesey, 2023) and were added using the rphylopic R package (Gearty & Jones, 2023). The silhouette was contributed as follows: “Elephas maximus” by T. Michael Keesey (2011; CC BY 3.0).

Figure 3. Elasticities of \(R_0\) with respect to matrix entries. Elasticities quantify the proportional effect on \(R_0\) of proportional changes in each vital rate. Panel (a) shows elasticities with respect to fertility rates (first row), and panel (b) shows elasticities with respect to survival probabilities (subdiagonal). Bars are grouped by vital rate, with colors indicating the five different models (original coarse model for the Nagarhole National Park population
and four 5-year aggregations of the model for the Periyar Reserve population). Elasticities are normalized to sum to 1 within each matrix, so the panels show how the influence on \(R_0\) is distributed across life-history components. The models in the figure legend are as described in caption of Table 1.


#Asian elephant
#"Elephas maximus" by  T. Michael Keesey (2011; CC BY 3.0).

# Not run: requires internet access

# uuid <- "7c9ab182-175d-4f02-96d0-09c1e5212bff"
# img <- get_phylopic(uuid = uuid, format = "raster", height = 691)
# Write image 
# png::writePNG(unclass(img), "inst/extdata/elephant.png")


img_file <- system.file("extdata", "elephant.png", package = "mpmaggregate")

if (!nzchar(img_file)) {
  img_file <- "../inst/extdata/elephant.png"
}

img <- png::readPNG(img_file)

# Compute elasticity matrices for each model (same dimensions as A)
# These are elasticities of R0 wrt to matrix entries
# normalized to sum to 1.
E_models <- lapply(models, function(A) {
  n <- nrow(A)
  matU <- A
  matU[1,] <- 0
  matF <- matrix(0,n,n)
  matF[1,] <- A[1,]
  E <- mpm_elasticity(matA=NULL, matF = matF, matU = matU, 
                  framework = "R0", normalize = TRUE)$elasticity
  as.matrix(E)
})

# ---- Helpers to extract Leslie elasticities from an elasticity matrix ----
extract_fertility_elast <- function(E) {
  as.numeric(E[1, ])
}

extract_survival_elast <- function(E) {
  n <- nrow(E)
  e_sub  <- E[cbind(2:n, 1:(n - 1))]  # elasticity of A[i+1, i]
  as.numeric(e_sub)
}

model_names <- names(models)

# Same colors as the vital-rate plot for consistency
cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e")

# Build elasticity-by-model matrices so:
# rows = models (5 bars), cols = vital rates (groupings)
F_el_by_model <- t(sapply(E_models, extract_fertility_elast))  # 5 x n
S_el_by_model <- t(sapply(E_models, extract_survival_elast))   # 5 x n-1

# Name columns (vital rates)
n <- nrow(models[[1]])
colnames(F_el_by_model) <- paste0("F", 1:n)
colnames(S_el_by_model) <- c(paste0("S", 1:(n - 1)))

# Keep only elasticities that are nonzero in at least one model
keep_F <- colSums(abs(F_el_by_model) > 0) > 0
keep_S <- colSums(abs(S_el_by_model) > 0) > 0

F_el_plot2 <- F_el_by_model[, keep_F, drop = FALSE]
S_el_plot2 <- S_el_by_model[, keep_S, drop = FALSE]

# Leave extra space on the right for legends that won't overlap bars
op <- par(
  mfrow = c(2, 1),
  mar = c(7, 4, 3, 7),
  oma = c(0, 0, 0, 0)
)

op <- par(
  mfrow = c(2, 1),
  mar = c(5, 4, 1.5, 7),
  oma = c(0, 0, 0, 0)
)

# Panel (a): Elasticities with respect to fertility rates
barplot(
  F_el_plot2,
  beside = TRUE,
  ylim = c(0, 0.20),
  col = cols,
  border = "white",
  las = 2,
  ylab = expression("Normalized elasticity of " * R[0] * ""),
  main = "Elasticities with respect to fertility rates"
)
mtext("(a)",
      side = 3,
      adj = 0,
      line = 0.2)
lims <-  par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x =  NULL,
  y = mydiff * .75,
  height =  .3 * mydiff * dim(img)[1] / 1024
)
par(xpd = NA)
legend(
  "topright",
  inset = c(-0.15, -0.05),
  legend = model_names,
  fill = cols,
  bty = "n",
  cex = 0.9
)
par(xpd = FALSE)

# Panel (b): Elasticities with respect to survival probabilities
barplot(
  S_el_plot2,
  beside = TRUE,
  ylim = c(0, .20),
  col = cols,
  border = "white",
  las = 2,
  ylab = expression("Normalized elasticity of " * R[0] * ""),
  main = "Elasticities with respect to survival probabilities"
)
mtext("(b)",
      side = 3,
      adj = 0,
      line = 0.2)
lims <-  par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x =  NULL,
  y = mydiff * .75,
  height =  .3 * mydiff * dim(img)[1] / 1024
)

par(xpd = NA)
legend(
  "topright",
  inset = c(-0.15, -0.05),
  legend = model_names,
  fill = cols,
  bty = "n",
  cex = 0.9
)

par(xpd = FALSE)

par(op)

How to read Figure 3. A larger elasticity (of \(R_0\) with respect to matrix entries) for a given vital rate means that small proportional changes to that vital rate would have a comparatively larger effect on net reproductive rate. For Asian elephants, elasticities are concentrated in survival probabilities rather than fertility rates, reflecting the strong influence of survivorship on \(R_0\). Aggregation methods yield similar elasticities. Late-life fertility rates and survival probabilities have higher elasticities in the Nagarhole National Park population than in the Periyar Reserve population.The silhouette is from PhyloPic (https://www.phylopic.org; Keesey, 2023) and were added using the rphylopic R package (Gearty & Jones, 2023). The silhouette was contributed as follows: “Elephas maximus” by T. Michael Keesey (2011; CC BY 3.0).

Comparing elasticities of \(\lambda\) and \(R_0\)

We compare elasticities of \(\lambda\) and \(R_0\) primarily through visual inspection, focusing on the relative distribution of elasticity mass between fertility and survival components of the matrix.

Across models, elasticity patterns under the \(\lambda\) and \(R_0\) frameworks are similar. In both cases, elasticities are dominated by survival transitions, consistent with the life history of long-lived species such as Asian elephants. The allocation of elasticity mass to fertility terms differs only slightly between the two frameworks.

Specifically, the total elasticity mass associated with fertility rates is 0.17 under the \(\lambda\) framework and 0.14 under the \(R_0\) framework. This small difference indicates that elasticities of \(\lambda\) place slightly greater weight on fertility than do elasticities of \(R_0\), although survival probability remains the dominant contributor in both cases.

#compare elasticities with respect to fertility rates between frameworks 
#fertility rates in lambda framework 
  sum(F_el_plot)/5
#> [1] 0.1742462
#fertility rates in R0 framework
  sum(F_el_plot2)/5
#> [1] 0.1439543

Summary

This vignette demonstrates how age-structured matrix population models (MPMs) defined with different projection intervals can be compared by aggregating the coarser model with leslie_aggregate(). Aggregating the annual Leslie MPM for the Asian elephant Periyar Reserve population to a 5-year Leslie MPM allowed us to directly compare it to a 5-year Leslie MPM of the Nagarhole National Park population.

Estimates of \(\lambda\) and \(R_0\) for the Periyar Reserve population indicate faster population growth and shorter generation time relative to the Nagarhole National Park population (Table 1). The lower population growth rate of the Nagarhole National Park population is attributable to lower fertility rates, as demonstrated by Figure 1. Elasticities indicate the survival probability dominates fertility rates as the most influential group of vital rates controlling \(\lambda\) and \(R_0\).

References

Chelliah, K., Bukka, H., and Sukumar, R. (2013). Modeling harvest rates and numbers from age and sex ratios: a demonstration for elephant populations. Biological Conservation, 165, 54–61. https://doi.org/10.1016/j.biocon.2013.05.008

Gearty, W., and Jones, L. A. (2023). rphylopic: An R package for fetching, transforming, and visualising PhyloPic silhouettes. Methods in Ecology and Evolution, 14, 2700–2708. https://doi.org/10.1111/2041-210X.14221

Goswami, V. R., Vasudev, D., and Oli, M. K. (2014). The importance of conflict-induced mortality for conservation planning in areas of human–elephant co-occurrence. Biological Conservation, 176, 191–198. https://doi.org/10.1016/j.biocon.2014.05.026

Hinrichsen, R. A. (2023). Aggregation of Leslie matrix models with application to ten diverse animal species. Population Ecology, 65(3), 146–166. https://doi.org/10.1002/1438-390X.12149

Hinrichsen, R. A., Yokomizo, H., and Salguero-Gómez, R. (2026). From theory to application: Elasticity-consistent aggregation of Leslie matrix population models for comparative demography. bioRxiv, preprint. https://doi.org/10.64898/2026.02.04.703802

Jones, O. R., Barks, P., Stott, I., James, T. D., Levin, S., Petry, W. K., Capdevila, P., Che-Castaldo, J., Jackson, J., Römer, G., Schuette, C., Thomas, C. C., and Salguero-Gómez, R. (2022). Rcompadre and Rage—Two R packages to facilitate the use of the COMPADRE and COMADRE databases and calculation of life-history traits from matrix population models. Methods in Ecology and Evolution, 13(12), 2700–2708. https://doi.org/10.1111/2041-210X.13792

Salguero-Gómez, R., Jones, O. R., Archer, C. R., Bein, C., de Buhr, H., Farack, C., Gottschalk, F., Hartmann, A., Henning, A., Hoppe, G., Römer, G., Ruoff, T., Sommer, V., Wille, J., Voigt, J., Zeh, S., Vieregg, D., Buckley, Y. M., Che-Castaldo, J., Hodgson, D., Scheuerlein, A., Caswell, H., and Vaupel, J. W. (2016). COMADRE: a global database of animal demography. Journal of Animal Ecology, 85, 371–384. https://doi.org/10.1111/1365-2656.12482