---
title: "Aggregating Asian Elephant Matrix Population Models for Comparative Demography with mpmaggregate"
author:
  - name: "Richard A. Hinrichsen"
    orcid: "0000-0003-0761-3005"
output:
  rmarkdown::html_vignette:
    number_sections: false
    mathjax: default
vignette: >
  %\VignetteIndexEntry{Aggregating Asian Elephant Matrix Population Models for Comparative Demography with mpmaggregate}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
  
**ORCID:** <https://orcid.org/0000-0003-0761-3005>
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```


## 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):

- Population/model ID: **249273** — a coarser model with a 5-year projection interval, 
Nagarhole National Park population, India (Chelliah et al. 2013)  
- Population/model ID: **249274** — a fine-grained model with a 1-year projection interval, 
Periyar Reserve population, India (Goswami et al. 2014)

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.

```{r packages}
library(mpmaggregate)
library(kableExtra)
library(rphylopic)
```
## Data acquisition (COMADRE via Rcompadre)

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

```{r fetch-comadre}

# 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
```{r leslie-aggregation}

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.

```{r leslie-metrics}

# 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

```{r define-models}
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\)

```{r lambda - r0 - table}

#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
  )
```


**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

- Aggregation under the \(\lambda\) framework prioritizes preservation of long-term population growth.
- Aggregation under the \(R_0\) framework emphasizes lifetime reproductive output.

## 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.


```{r vital-rates-barplot, fig.height=8, fig.width=10}

#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.


```{r lambda-elasticities-barplot, fig.height=8, fig.width=10}
#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.


```{r R0-elasticities-barplot, fig.height=8, fig.width=10}

#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.

```{r lambda-vs.-R0-elasticities}
#compare elasticities with respect to fertility rates between frameworks 
#fertility rates in lambda framework 
  sum(F_el_plot)/5
#fertility rates in R0 framework
  sum(F_el_plot2)/5
```

## 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
