## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----packages-----------------------------------------------------------------
library(mpmaggregate)
library(kableExtra)
library(rphylopic)

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

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

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


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

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

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

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

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

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