--- title: "Aggregating a Leslie-to-Leslie MPM" author: - name: "Richard A. Hinrichsen" orcid: "0000-0003-0761-3005" output: rmarkdown::html_vignette: number_sections: false mathjax: default vignette: > %\VignetteIndexEntry{Aggregating a Leslie-to-Leslie MPM} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- **ORCID:** ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) ``` ## Overview Matrix population models (MPMs) are constructed at different age or stage resolutions, which complicates demographic comparisons across populations. For age-structured populations, such complications arise when a Leslie matrix model (MPM1) with a small age-class width is compared with a coarser Leslie matrix model (MPM2). In such a case, aggregation can be used to group fine age classes of MPM1 into a smaller number of wider age classes that are compatible with MPM2, thus allowing direct comparison of the demographic features. Aggregation must proceed carefully because certain demographic parameters are known to change with aggregation. In this vignette, we demonstrate Leslie-to-Leslie MPM aggregation using the `mpmaggregate` package. Leslie-to-Leslie MPM aggregation proceeds by reducing the number of age classes of the Leslie MPM while retaining its Leslie form. Unlike general-to-general MPM aggregation, Leslie-to-Leslie MPM aggregation changes the projection interval to ensure that it is equal to the age-class width. As a worked example, we use an age-structured matrix population model for the yellow-bellied marmot *Marmota flaviventris* obtained from the COMADRE database (Jones et al. 2022), which archives demographic data originally published by Ozgul et al. (2009). We focus on aggregation using `leslie_aggregate()`, and we compare four combinations of demographic framework and aggregation criterion, highlighting what each approach preserves and how the aggregated models differ in interpretation. ## 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 `expm` is used to raise a matrix to a power, `knitr` is used for report generation, and `kableExtra` is used for table formatting. ```{r packages} library(mpmaggregate) library(kableExtra) library(rphylopic) library(expm) ``` ## Data acquisition To demonstrate Leslie-to-Leslie MPM aggregation, we use an age-structured MPM for *Marmota flaviventris* (MatrixID **249455**) from the COMADRE database (Ozgul et al. 2009). We aggregate this Leslie MPM with 10 age classes to a Leslie MPM with 5 age classes. Age classes 1-2 are combined to form a new age class 1, age classes 3-4 are combined to form a new age class 2, etc.. The yellow-bellied marmot *Marmota flaviventris* is a large, stout-bodied ground squirrel (rodent) native to the mountainous regions of western North America. ```{r} # Leslie population projection matrix for Marmota flaviventris # MatrixID = 249455 # Not run: requires Rcompadre and internet access # Matrices can be retrieved with: # library(Rcompadre) # comadre <- cdb_fetch("comadre",version = "4.23.3.1") # mpm <- comadre[comadre$MatrixID == 249455, ] # matA <- matA(mpm)[[1]] # matU <- matU(mpm)[[1]] # matF <- matF(mpm)[[1]] # Matrices are defined locally to avoid requiring the use of the internet and Rcompadre # Leslie population projection matrix matA <- matrix(c( 0.000000, 0.258040, 0.627254, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.545125, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.498625, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.592875, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000 ), nrow = 10, byrow = TRUE) # Reproduction matrix matF <- matrix(c( 0, 0.25804, 0.627254, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ), nrow = 10, byrow = TRUE) # Survival/transition matrix matU <- matrix(c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.545125,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.498625,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.592875,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6435, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6435, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6435, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6435, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6435, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6435, 0 ), nrow = 10, byrow = TRUE) # Sanity check: A = U + F (after folding any C into F) stopifnot(all.equal(unname(matA), unname(matU + matF))) ``` ## Helper functions for demographic quantities We compute: - \( \lambda \): spectral radius of \( \mathbf{A} \) - \( R_0 \): spectral radius of \( \mathbf{F} \left( \mathbf{I} - \mathbf{U} \right)^{-1} \) - \( T_a \): generation time (years) under the \( \lambda \) framework - \( T_c \): cohort generation time (years) under the \( R_0 \) framework ```{r helpers} `%||%` <- function(x, y) if (!is.null(x)) x else y # Return R0, the net reproductive rate get_R0 <- function(U, F) { I <- diag(nrow(U)) N <- solve(I - U) K <- F %*% N spectral_radius(K) } # Return generation time get_Ta <- function(U, F) { generation_time(F, U, framework = "lambda")$generation_time } # Return cohort generation time get_Tc <- function(U, F) { generation_time(F, U, framework = "R0")$generation_time } ``` ## Aggregate in four ways We aggregate the model using the following four combinations of framework and criterion: 1. `framework = "lambda"`, `criterion = "standard"` 2. `framework = "lambda"`, `criterion = "elasticity"` 3. `framework = "R0"`, `criterion = "standard"` 4. `framework = "R0"`, `criterion = "elasticity"` In the \(\lambda\) framework, standard aggregation follows the method of Hooley (2000), while elasticity-consistent aggregation follows Bienvenu et al. (2017). These methods are extended to Leslie-to-Leslie aggregation (Hinrichsen, 2023; Hinrichsen et al., 2026). Both approaches maintain consistency of the asymptotic growth rate \(\lambda\) and the stable age distribution. Elasticity-consistent aggregation additionally maintains consistency of the reproductive values and the elasticities of \(\lambda^k\) with respect to entries of \( \mathbf{A}^k \), where \(k = n/m\). We extend these Leslie-to-Leslie approaches to the \(R_0\) framework, where standard aggregation maintains the consistency of the net reproductive rate \(R_0\) and cohort stable age distribution. Elasticity-consistent aggregation in this framework additionally maintains the consistency of the cohort reproductive values and elasticities of \(R_0\) with respect to matrix entries of \( \mathbf{A}_k \), which is a k-step-ahead version of \( \mathbf{A}\) in the \(R_0\) framework. ```{r aggregate} # Aggregate into 5 x 5 matrices (m = 5) ngroups <- 5 agg1 <- leslie_aggregate( matA = matA, ngroups = ngroups, framework = "lambda", criterion = "standard" ) agg2 <- leslie_aggregate( matA = matA, ngroups = ngroups, framework = "lambda", criterion = "elasticity" ) agg3 <- leslie_aggregate( matA = matA, ngroups = ngroups, framework = "R0", criterion = "standard" ) agg4 <- leslie_aggregate( matA = matA, ngroups = ngroups, framework = "R0", criterion = "elasticity" ) # Split a Leslie matrix into U and F components (Leslie form: fertility in row 1) split_matA <- function(A) { n <- nrow(A) U <- A U[1, ] <- 0 F <- matrix(0, n, n) F[1, ] <- A[1, ] list(U = U, F = F, A = A) } m0 <- list(U = matU, F = matF, A = matA) m1 <- split_matA(agg1$matAk_agg) m2 <- split_matA(agg2$matAk_agg) m3 <- split_matA(agg3$matAk_agg) m4 <- split_matA(agg4$matAk_agg) eff <- c( NA, agg1$effectiveness, agg2$effectiveness, agg3$effectiveness, agg4$effectiveness ) models <- list( Original = m0, Agg_lambda_std = m1, Agg_lambda_elast = m2, Agg_R0_std = m3, Agg_R0_elast = m4 ) model_labels <- c( Original = "Original", Agg_lambda_std = "Agg λ + standard", Agg_lambda_elast = "Agg λ + elasticity", Agg_R0_std = "Agg R0 + standard", Agg_R0_elast = "Agg R0 + elasticity" ) # Basic checks stopifnot(all(sapply(models, function(m) !is.null(m$A)))) # Show aggregated matrices print("Agg λ + standard"); m1$A print("Agg λ + elasticity"); m2$A print("Agg R0 + standard"); m3$A print("Agg R0 + elasticity");m4$A ``` **Important.** Unlike general-to-general MPM aggregation, Leslie-to-Leslie MPM aggregation cannot produce survival probabilities greater than 1. Furthermore, within each framework, the survival probabilities are the same under standard and elasticity-consistent aggregation, while the fertility schedules can differ between the standard and elasticity-consistent aggregated Leslie matrices (Hinrichsen et al., 2026). ## Compare \(\lambda\), \(R_0\), \(T_a\), and \(T_c\) in a table ```{r results-table alt} # Adjust lambdas so that they are all per capita growth rate over a year # Adjust generation times so they all use units of 1 year adj = c(1,2,2,2,2) lambda = sapply(models, function(m) spectral_radius(m$A)) lambda = lambda^(1/adj) R0 = sapply(models, function(m) get_R0(m$U, m$F)) Ta = sapply(models, function(m) get_Ta(m$U, m$F))*adj Tc = sapply(models, function(m) get_Tc(m$U, m$F))*adj results <- data.frame( Model = unname(model_labels[names(models)]), lambda = lambda, R0 = R0, Ta = Ta, Tc = Tc, Effectiveness = eff, row.names = NULL, stringsAsFactors = FALSE ) n_rows <- nrow(results) knitr::kable( results, col.names = c( "Model", "λ", "R0", "Ta (years)", "Tc (years)", "ρ2" ), digits = 3, caption = paste0( "Table 1. Comparison of demographic properties for the ", "yellow-bellied marmot Marmota flaviventris (COMADRE MatrixID 249455) ", "using the original 10 × 10 Leslie matrix and four aggregated Leslie matrices ", "collapsed to 5 age groups. 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." ), 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( "", "Note. Population growth rates (λ) are all transformed to ", "represent growth over 1 year. Generation times have all been transformed to ", "represent generation time measured in years." ), footnote_as_chunk = TRUE, escape = FALSE ) ``` **How to read Table 1.** The first row reports demographic quantities computed from the original Leslie matrix for *Marmota flaviventris* (Ozgul et al. 2009), while the remaining rows report estimates obtained from aggregated Leslie matrices. ### Interpretation notes - Aggregation under the **\(\lambda\) framework** prioritizes consistency of stable population growth rate. - Aggregation under the **\(R_0\) framework** prioritizes consistency of net reproductive rate. ## Elasticities of \(\lambda\): making them comparable across models To make the elasticities of the original model comparable with those of the aggregated models, we map the original elasticity matrix to the aggregated age structure as \[ \mathbf{E}^{\mathrm{true}}_{\mathrm{agg}} = \mathbf{P}\, \mathbf{E}\!\bigl(\mathbf{A}^{k}\bigr)\, \mathbf{P}^{\top}, \] where \( \mathbf{E}\!\bigl(\mathbf{A}^{k}\bigr) \) is the elasticity of \( \lambda^{k} \) with respect to the entries of the original \( \mathbf{A}^{k} \), and \( \mathbf{P} \) is the partitioning matrix induced by `groups`. Elasticities for each aggregated model are then computed directly from its aggregated population projection matrix, and the results for all models are plotted together. ```{r elasticities-lambda^k} # Age groups used for aggregation in example groups <- list(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10)) # Partitioning matrix induced by the age groups. # This matrix is constructed internally by leslie_aggregate(), # but we construct it explicitly here so it can be used to # map quantities from the original age-class space to the # aggregated space. P <- mpm_partition(groups = groups, n = nrow(matA)) # Dimensions of the original and aggregated matrices m <- ngroups n <- nrow(matA) k <- n / m matAk <- matA %^% k # Elasticity of lambda^k with respect to entries of A^k E_A <- mpm_elasticity(matA = matAk, framework = "lambda")$elasticity E_list <- list( "Original" = P %*% E_A %*% t(P), "Agg λ + standard" = mpm_elasticity(matA = m1$A, framework = "lambda")$elasticity, "Agg λ + elasticity" = mpm_elasticity(matA = m2$A, framework = "lambda")$elasticity, "Agg R0 + standard" = mpm_elasticity(matA = m3$A, framework = "lambda")$elasticity, "Agg R0 + elasticity" = mpm_elasticity(matA = m4$A, framework = "lambda")$elasticity ) to_long <- function(M, model_name) { idx <- which(M != 0, arr.ind = TRUE) data.frame( model = model_name, row = idx[, 1], col = idx[, 2], entry = paste(idx[, 1], idx[, 2], sep = ","), elasticity = M[idx], stringsAsFactors = FALSE ) } elast_df <- do.call(rbind, Map(to_long, E_list, names(E_list))) elast_df <- elast_df[elast_df$elasticity > 0, ] elast_df <- elast_df[order(elast_df$row, elast_df$col), ] models_order <- names(E_list) entries <- unique(elast_df$entry) elast_mat <- sapply(models_order, function(m) elast_df$elasticity[elast_df$model == m]) elast_mat <- t(elast_mat) rownames(elast_mat) <- models_order colnames(elast_mat) <- entries ``` ### Plot: elasticities of \( \lambda^k \) by matrix entry Each x-axis category corresponds to a matrix entry \([i,j]\) that is nonzero in at least one of the five \(A\) matrices. Within each category, we plot elasticities for the original (“true aggregated”) model and the four aggregated models. **Figure 1. Elasticities of \( \lambda^k \) with respect to matrix entries.** Elasticities quantify the proportional change in \( \lambda^k \) resulting from proportional changes in entries of the entries of \( \mathbf{A}^{k} \) . Bars are grouped by matrix entry, with colors indicating the five models. Because elasticities sum to 1 within each matrix, the figure illustrates how the influence on \( \lambda^k \) is distributed across life-history components. The models shown in the legend correspond to those described in Table 1. ```{r elasticity-plot, fig.width=12, fig.height=6} # Not run: requires internet access # Marmot 2012 Feb 2 by T. Michael Keesey, CC0 1.0, Marmota monax (Linnaeus 1758) # uuid <- "eee50efb-40dc-47d0-b2cb-52a14a5e0e51" # img <- get_phylopic(uuid = uuid, format = "raster", height = 869) # Write image # png::writePNG(unclass(img), "inst/extdata/marmot.png") # Load image locally img_file <- system.file("extdata", "marmot.png", package = "mpmaggregate") if (!nzchar(img_file)) { img_file <- "../inst/extdata/marmot.png" } img <- png::readPNG(img_file) rc <- do.call(rbind, strsplit(colnames(elast_mat), ",")) r_vec <- as.integer(rc[, 1]) c_vec <- as.integer(rc[, 2]) make_entry_label <- function(r, c) { if (r == 1) { bquote(F[1, .(c)]) } else { bquote(U[.(r), .(c)]) } } entry_labels <- mapply(make_entry_label, r_vec, c_vec, SIMPLIFY = FALSE) cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e") op <- par(mar = c(7, 6, 4, 2)) bp <- barplot( elast_mat, beside = TRUE, ylim = c(0, .35), col = cols, border = NA, ylab = expression("Elasticity of " * lambda^k), xaxt = "n", space = c(0.2, 1) ) axis(side = 1, at = colMeans(bp), labels = entry_labels, las = 2, cex.axis = 0.9) legend("topleft", legend = rownames(elast_mat), fill = cols, bty = "n", inset = 0.02) lims <- par("usr") mydiff <- lims[4] - lims[3] add_phylopic_base( img = img, x = NULL, y = mydiff * .9, height = .2 * mydiff * dim(img)[1] / 1024 ) par(op) ``` **How to read Figure 1.** Each x-axis label corresponds to a vital rate represented by a matrix entry in the aggregated space. Fertility rate, \(F[1,i]\) denotes reproduction of age class \(i\) over the projection interval. Survival probability, \(U[i+1,i]\) gives the probability of surviving and transitioning from age class \(i\) to age class \(i+1\) over the projection interval. There is exact agreement between the elasticities of models “Original” and “Agg \(\lambda\) + elasticity,” which is by design. 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: *Marmota monax* by T. Michael Keesey (2012; CC0 1.0). ## Elasticities of \(R_0\): making them comparable across models To make the elasticities of the original model comparable with those of the aggregated models, we map the original elasticity matrix to the aggregated age structure as \[ \mathbf{E}^{\mathrm{true}}_{\mathrm{agg}} = \mathbf{P}\, \mathbf{E}\!\left(\mathbf{A}_{k}\right)\, \mathbf{P}^{\top}, \] where \( \mathbf{E}\!\left(\mathbf{A}_{k}\right) \) is the normalized elasticity of \( R_0 \) with respect to the entries of the \( k \)-step-ahead population projection matrix \( \mathbf{A}_{k} \) in the \( R_0 \) framework, and \( \mathbf{P} \) is the partitioning matrix induced by `groups`. Elasticities for each aggregated model are then computed directly from its corresponding \( \mathbf{U}_{k} \) and \( \mathbf{F}_{k} \), and the results for all models are plotted together. ```{r elasticities-r0} # Age groupings used for aggregation (reintroduced here so this section # can be read independently) groups <- list(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10)) # Partitioning matrix induced by the age groups. # Constructed internally by leslie_aggregate(), but reproduced here # because it is used to map quantities to aggregated space. P <- mpm_partition(groups = groups, n = nrow(matA)) # Dimensions of the original and aggregated matrices # Reintroduced here for clarity m <- ngroups n <- nrow(matA) k <- n / m R0 <- spectral_radius(matF %*% solve(diag(n) - matU)) # Reference operator for aggregation eigenstructure in the R0 framework # Divide by R0 so that it has Leslie matrix form A_R0ref <- (matF + R0 * matU) / R0 Uk <- expm::`%^%`(matU, k) Rk <- (expm::`%^%`(A_R0ref, k) - Uk) * R0 Ak <- Uk + Rk # Elasticity of R0 with respect to entries of Ak (closest to what aggregated matrices yield) E_A <- mpm_elasticity(matF = Rk, matU = Uk, framework = "R0", normalize = TRUE)$elasticity E_list2 <- list( "Original" = P %*% E_A %*% t(P), "Agg λ + standard" = mpm_elasticity(matF = m1$F, matU = m1$U, framework = "R0", normalize = TRUE)$elasticity, "Agg λ + elasticity" = mpm_elasticity(matF = m2$F, matU = m2$U, framework = "R0", normalize = TRUE)$elasticity, "Agg R0 + standard" = mpm_elasticity(matF = m3$F, matU = m3$U, framework = "R0", normalize = TRUE)$elasticity, "Agg R0 + elasticity" = mpm_elasticity(matF = m4$F, matU = m4$U, framework = "R0", normalize = TRUE)$elasticity ) elast_df2 <- do.call(rbind, Map(to_long, E_list2, names(E_list2))) elast_df2 <- elast_df2[elast_df2$elasticity > 0, ] elast_df2 <- elast_df2[order(elast_df2$row, elast_df2$col), ] models_order2 <- names(E_list2) entries2 <- unique(elast_df2$entry) elast_mat2 <- sapply(models_order2, function(m) elast_df2$elasticity[elast_df2$model == m]) elast_mat2 <- t(elast_mat2) rownames(elast_mat2) <- models_order2 colnames(elast_mat2) <- entries2 ``` ### Plot: elasticities of \(R_0\) by matrix entry Each x-axis category corresponds to a matrix entry \([i,j]\) that is nonzero in at least one of the five matrices. Within each category, we plot elasticities for the original (“true aggregated”) model and the four aggregated models. **Figure 2. Elasticities of \(R_0\) with respect to matrix entries.** Elasticities quantify the proportional change in \(R_0\) resulting from proportional changes in individual vital rates. Bars are grouped by matrix entry, with colors indicating the five models. The elasticities of \(R_0\) are normalized to sum to 1, so each elasticity can be interpreted as the fraction of the total contribution to \(R_0\) attributable to that matrix entry. ```{r elasticity-r0-plot, fig.width=12, fig.height=6} # Not run: requires internet access # Marmot 2012 Feb 2 by T. Michael Keesey, CC0 1.0, Marmota monax (Linnaeus 1758) # uuid <- "eee50efb-40dc-47d0-b2cb-52a14a5e0e51" # img <- get_phylopic(uuid = uuid, format = "raster", height = 869) # Write image # png::writePNG(unclass(img), "inst/extdata/marmot.png") # Load the image file locally img_file <- system.file("extdata", "marmot.png", package = "mpmaggregate") #For development if (!nzchar(img_file)) { img_file <- "../inst/extdata/marmot.png" } img <- png::readPNG(img_file) rc <- do.call(rbind, strsplit(colnames(elast_mat2), ",")) r_vec <- as.integer(rc[, 1]) c_vec <- as.integer(rc[, 2]) entry_labels <- mapply(make_entry_label, r_vec, c_vec, SIMPLIFY = FALSE) cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e") op <- par(mar = c(7, 6, 4, 2)) bp <- barplot( elast_mat2, beside = TRUE, ylim = c(0, .35), col = cols, border = NA, ylab = expression("Normalized elasticity of " * R[0]), xaxt = "n", space = c(0.2, 1) ) axis(side = 1, at = colMeans(bp), labels = entry_labels, las = 2, cex.axis = 0.9) legend("topleft", legend = rownames(elast_mat2), fill = cols, bty = "n", inset = 0.02) lims <- par("usr") mydiff <- lims[4] - lims[3] add_phylopic_base( img = img, x = NULL, y = mydiff * .9, height = .2 * mydiff * dim(img)[1] / 1024 ) par(op) ``` **How to read Figure 2.** Each x-axis label corresponds to a vital rate represented by a matrix entry in the aggregated space. Fertility rate, \(F[1,i]\) denotes reproduction of age class \(i\) over the projection interval. Survival probability, \(U[i,i+1]\) gives the probability of surviving and transitioning from age class \(i\) to age class \(i+1\) over the projection interval. There is exact agreement between the elasticities of models “Original” and “Agg \(R_0\) + elasticity,” which is by design. The silhouette is from PhyloPic (https://www.phylopic.org; Keesey, 2023) and wereadded using the rphylopic R package (Gearty & Jones, 2023). The silhouette was contributed as follows: *Marmota monax* by T. Michael Keesey (2012; CC0 1.0). ### Comparing elasticities of \(\lambda^k\) and \(R_0\) Within each framework, true elasticities and those derived from aggregated matrices agree qualitatively, regardless of aggregation method. Both frameworks identify early survival and early fertility as dominant contributors \(\lambda^k\) and \(R_0\). Overall, fertility rates play a larger role in determining \(R_0\) than \(\lambda^k\), accounting for approximately 51\% versus 45\% of total elasticity mass, respectively. Across analyses, early survival probability (the transition \(U[2,1]\)) is the single most influential vital rate. ```{r fertility-elasticity-mass} # Mass for elasticity of R0 w.r.t. fertilities sum(elast_mat2[, 1:5]) / 5 # Mass for elasticity of lambda^k w.r.t. fertilities sum(elast_mat[, 1:5]) / 5 ``` ## References Bienvenu, F., Akçay, E., Legendre, S., and McCandlish, D. M. (2017). The genealogical decomposition of a matrix population model with applications to the aggregation of stages. *Theoretical Population Biology*, 115, 69–80. https://doi.org/10.1016/j.tpb.2017.04.002 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 Hooley, D. E. (2000). Collapsed matrices with (almost) the same eigenstuff. *The College Mathematics Journal*, 31(4), 297–299. https://doi.org/10.1080/07468342.2000.11974162 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(4), 770–781. https://doi.org/10.1111/2041-210X.13792 Ozgul, A., Oli, M. K., Armitage, K. B., Blumstein, D. T., and Van Vuren, D. H. (2009). Influence of local demography on asymptotic and transient dynamics of a yellow-bellied marmot metapopulation. *The American Naturalist*, 173(4), 517–530. https://doi.org/10.1086/597225