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