## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----projection matrix, echo=FALSE, out.width="50%", fig.cap="  "-------------
knitr::include_graphics("figures/Slide2.PNG")

## ----lambda-r0-figure, echo=FALSE, out.width="90%", fig.cap="Figure 1. Interpretation of λ and R0 frameworks"----
knitr::include_graphics("figures/Slide1.PNG")

## -----------------------------------------------------------------------------
library(mpmaggregate)
library(expm)
library(knitr)
library(kableExtra)

## ----table-preserved, echo=FALSE, results='asis'------------------------------

preserved_tbl <- data.frame(
  Framework = c("lambda", "lambda", "R0", "R0"),
  Criterion = c("standard", "elasticity-consistent", "standard", "elasticity-consistent"),
  Preserved = c(
    "<em>λ</em> (stable growth rate)<br>
     Stable stage distribution<br>
     Interstage flow matrix",

    "<em>λ</em> (stable growth rate)<br>
     Stable stage distribution<br>
     Reproductive values<br>
     Elasticities of <em>λ</em><br>
     Generation time (<em>T</em><sub>a</sub>)",

    "<em>R</em><sub>0</sub> (net reproductive rate)<br>
     Cohort stable stage distribution<br>
     Cohort interstage flow matrix",

    "<em>R</em><sub>0</sub> (net reproductive rate)<br>
     Cohort stable stage distribution<br>
     Cohort reproductive values<br>
     Elasticities of <em>R</em><sub>0</sub><br>
     Cohort generation time (<em>T</em><sub>c</sub>)"
  ),
  check.names = FALSE
)

knitr::kable(
  preserved_tbl,
  format = "html",
  escape = FALSE,
  align = c("l", "l", "l"),
  col.names = c("Framework", "Criterion", "Consistent Demographic Quantities"),
caption = paste0(
  "<strong>Table 1.</strong> Demographic quantities that remain consistent under ",
  "<code>mpm_aggregate()</code> in the <em>λ</em> and <em>R</em><sub>0</sub> frameworks ",
  "for both the standard and elasticity-consistent aggregation criteria. ",
  "Under the <em>λ</em> framework, aggregation maintains consistency in stable ",
  "growth properties, whereas under the <em>R</em><sub>0</sub> framework it maintains ",
  "consistency in cohort-based reproductive quantities. ",
  "Interstage flow weights the columns of <strong>A</strong> by their respective ",
  "stable stage densities (Yokomizo et al. 2024). ",
  "Here, “cohort” denotes the <em>R</em><sub>0</sub>-framework analogues of the ",
  "corresponding <em>λ</em>-framework demographic quantities. ",
  "Notice the symmetry between the <em>λ</em> and <em>R</em><sub>0</sub> frameworks."
)
) |>
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "condensed")
  ) |>
  column_spec(3, width = "7cm") |>
  row_spec(0, extra_css = "border-bottom: 2px solid black;") |>
  footnote(
    general_title = "",
    general = paste0(
      "<div style='border-top: 1px solid black; padding-top: 6px;'>",
      "<span style='font-size: 90%;'>",
      "<strong>Note.</strong> <code>mpm_aggregate()</code> retains the projection ",
      "interval of the original matrix population model.",
      "</span></div>"
    ),
    footnote_as_chunk = TRUE,
    escape = FALSE
  )

## -----------------------------------------------------------------------------
# Example matrices used in the general-to-general MPM aggregation examples below
#population projection matrix for Sclerocarya	birrea
#MatrixID = 247940

# Not run: requires Rcompadre and internet access

#these matrices can be retrieved with:
#library(Rcompadre)
#compadre <- cdb_fetch("compadre")
#mpm <- compadre[compadre$MatrixID == 247940, ]
#A <- matA(mpm)[[1]]
#U <- matU(mpm)[[1]]
#F <- matF(mpm)[[1]]


A <- matrix(
  c(0.0, 0.0000, 12.46,
    0.4, 0.9160,  0.00,
    0.0, 0.0371,  0.99),
  nrow = 3,
  byrow = TRUE
)
#survival/transitions matrix
U <- matrix(
  c(0.0, 0.0000,  0.00,
    0.4, 0.9160,  0.00,
    0.0, 0.0371,  0.99),
  nrow = 3,
  byrow = TRUE
)
#reproduction matrix
F <- matrix(
  c(0.0, 0.0000, 12.46,
    0.0, 0.0000,  0.00,
    0.0, 0.0000,  0.00),
  nrow = 3,
  byrow = TRUE
)

# Stage aggregation used in later examples:
# group vegetative + reproductive stages together,
# leaving the seedling stage separate
groups <- list(
  c(1),
  c(2,3)
)

## -----------------------------------------------------------------------------

# General-to-general MPM aggregation
# in the lambda framework, standard aggregation
res_lambda_std <- mpm_aggregate(
  matA = A,
  groups = groups,
  framework = "lambda",
  criterion = "standard"
)

# Aggregated projection matrix
res_lambda_std$matA_agg

# Effectiveness of aggregation
res_lambda_std$effectiveness

# ---- Check: lambda consistency ----

# Expected
spectral_radius(A)

# Lambda of aggregated model
spectral_radius(res_lambda_std$matA_agg)

# ---- Check: stable stage distribution consistency ----

# P is the partitioning matrix that groups stages
P <- mpm_partition(groups=groups,n=3)

# w is the stable stage distribution
w <- stable_stage(A=A)

# Expected
c(P%*%w)

# stable stage distribution of aggregated model
w_agg <-stable_stage(res_lambda_std$matA_agg)
w_agg

# ---- Check: interstage flows consistency ----

# Expected
P%*%A%*%diag(w)%*%t(P)

# Interstage flows of aggregated model
res_lambda_std$matA_agg%*%diag(w_agg)

## -----------------------------------------------------------------------------

# General-to-general MPM aggregation
# in the lambda framework, elasticity-consistent aggregation
res_lambda_el <- mpm_aggregate(
  matA = A,
  matF = F,
  matU = U,
  groups = groups,
  framework = "lambda",
  criterion = "elasticity"
)

# Aggregated matrix
res_lambda_el$matA_agg

# Effectiveness of aggregation
res_lambda_el$effectiveness

# ---- Check: lambda consistency ----

# Expected
spectral_radius(A)

# Lambda of aggregated model
spectral_radius(res_lambda_el$matA_agg)

# ---- Check: stable stage distribution consistency ----

# P is the partitioning matrix
P <- mpm_partition(groups=groups,n=3)

# w us the stable stage distribution 
w <- stable_stage(A=A)

# Expected
c(P%*%w)

# Stable stage distribution of aggregated model
stable_stage(res_lambda_el$matA_agg)

# ---- Check: reproductive values consistency ----

# v is reproductive values (dominant left eigenvector)
v = reproductive_values(A)

# Expected
c(solve(P%*%diag(w)%*%t(P))%*%P%*%(w*v))

# Reproductive values of aggregated model
reproductive_values(res_lambda_el$matA_agg)

# ---- Check: elasticities consistency ----

# Elasticity matrix
EA <- mpm_elasticity(matA=A,framework="lambda")$elasticity

# Expected
P%*%EA%*%t(P)

# Elasticities of aggregated model
mpm_elasticity(matA=res_lambda_el$matA_agg,framework="lambda")$elasticity

# ---- Check: generation time consistency ----

# Expected
generation_time(matF=F,matU=U,framework="lambda")$generation_time

# Generation time from aggregated model
generation_time(matF=res_lambda_el$matF_agg,
                matU=res_lambda_el$matU_agg,
                framework="lambda")$generation_time

## -----------------------------------------------------------------------------

# General-to-general MPM aggregation
# in the R0 framework, standard aggregation
res_R0_std <- mpm_aggregate(
  matU = U,
  matF = F,
  groups = groups,
  framework = "R0",
  criterion = "standard"
)

# Aggregated projection matrix
res_R0_std$matA_agg

# Effectiveness of aggregation
res_R0_std$effectiveness

# ---- Check: R0 consistency ----

# Expected
spectral_radius(F%*%solve(diag(3)-U))

# R0 from aggregated matrix
spectral_radius(res_R0_std$matF_agg%*%solve(diag(2)-res_R0_std$matU_agg))

# ---- Check: cohort stable stage distribution consistency ----

# w is the cohort stable stage distribution from a next generation matrix
w <- stable_stage(solve(diag(3)-U)%*%F)

# Expected for consistency
c(P%*%stable_stage(solve(diag(3)-U)%*%F))

# Aggregated cohort stable stage distribution
w_agg <- stable_stage(solve(diag(2)-res_R0_std$matU_agg)%*%res_R0_std$matF_agg)
w_agg

# ---- Check: cohort interstage flows  consistency ----

# Expected
P%*%A%*%diag(w)%*%t(P)

# Aggregated interstage flows
res_R0_std$matA_agg%*%diag(w_agg)

## -----------------------------------------------------------------------------

# General-to-general MPM aggregation
# in the R0 framework, elasticity-consistent aggregation
res_R0_el <- mpm_aggregate(
  matU = U,
  matF = F,
  groups = groups,
  framework = "R0",
  criterion = "elasticity"
)

# Aggregated projection matrix
res_R0_el$matA_agg

# Effectiveness of aggregation
res_R0_el$effectiveness

# ---- Check: R0 consistency ----

# Expected R0
R0 <- spectral_radius(F%*%solve(diag(3)-U))
R0

# R0 from the aggregated model
spectral_radius(res_R0_el$matF_agg%*%solve(diag(2)-res_R0_el$matU_agg))

# ---- Check: cohort stable stage distribution consistency ----

# Get partitioning matrix
P <- mpm_partition(groups=groups,n=3)

# Construct cohort projection matrix
Aref = F + U*R0

# Cohort stable stage distribution
w <- stable_stage(A=Aref)

# Expected for consistency
c(P%*%w)

# Aggregated cohort projection matrix
Aref_agg = res_R0_el$matF_agg+ R0*res_R0_el$matU_agg

# Aggregated cohort stable stage distribution
stable_stage(Aref_agg)

# ---- Check: cohort reproductive values consistency ----
v = reproductive_values(Aref)

# Expected
c(solve(P%*%diag(w)%*%t(P))%*%P%*%(w*v))

# Cohort reproductive values for aggregated matrix
reproductive_values(Aref_agg)

# ---- Check: elasticities of R0 consistency ----
# Here elasticities are normalized to sum to 1

# Original elasticities
EA <- mpm_elasticity(matF=F,matU=U,framework="R0")$elasticity

# Expected for consistency
P%*%EA%*%t(P)

# Elasticities of R0 from aggregated matrix
mpm_elasticity(matF=res_R0_el$matF_agg,matU=res_R0_el$matU_agg,framework="R0")$elasticity

# ---- Check: cohort generation time consistency ----

# Expected (from original model)
generation_time(matF=F,matU=U,framework="R0")$generation_time

# Cohort generation time from aggregated model
generation_time(matF = res_R0_el$matF_agg, 
                matU = res_R0_el$matU_agg,
                framework="R0")$generation_time

## ----table-leslie-preserved, echo=FALSE, results='asis'-----------------------
preserved_leslie_tbl <- data.frame(
  Framework = c("lambda", "lambda", "R0", "R0"),
  Criterion = c("standard", "elasticity-consistent", "standard", "elasticity-consistent"),
  Preserved = c(
    "<em>λ</em> (stable growth rate)<br>
     Stable age distribution<br>
     Interstage flow matrix of <strong>A</strong><sup><em>k</em></sup>",

    "<em>λ</em> (stable growth rate)<br>
     Stable age distribution<br>
     Reproductive values<br>
     Elasticities of <em>λ</em><sup><em>k</em></sup> with respect to the entries of <strong>A</strong><sup><em>k</em></sup><br>
     Generation time (<em>T</em><sub>a</sub>) of <strong>A</strong><sup><em>k</em></sup>",

    "<em>R</em><sub>0</sub> (net reproductive rate)<br>
     Cohort stable age distribution<br>
     consistent of <strong>A</strong><sub><em>k</em></sub><sup></sup>",

    "<em>R</em><sub>0</sub> (net reproductive rate)<br>
     Cohort stable age distribution<br>
     Cohort reproductive values<br>
     Elasticities of <em>R</em><sub>0</sub> with respect to <strong>A</strong><sub><em>k</em></sub><sup></sup><br>
     Cohort generation time (<em>T</em><sub>c</sub>) for <strong>A</strong><sub><em>k</em></sub><sup></sup>"
  ),
  check.names = FALSE
)

knitr::kable(
  preserved_leslie_tbl,
  format = "html",
  escape = FALSE,
  align = c("l", "l", "l"),
  col.names = c("Framework", "Criterion", "Consistent Demographic Quantities"),
caption = paste0(
  "<strong>Table 2.</strong> Summary of the demographic properties that remain consistent under ",
  "<code>leslie_aggregate()</code> when reducing a Leslie matrix from dimensionality ",
  "<em>n</em> to <em>m</em>. The quantities the remain consistent depend on both the demographic framework ",
  "(long-term growth under <em>λ</em> or lifetime reproductive output under ",
  "<em>R</em><sub>0</sub>) and the aggregation criterion (standard versus ",
  "elasticity-consistent). The integer <em>k</em> = <em>n</em>/<em>m</em> defines the ",
  "temporal rescaling by aggregating a Leslie model. Quantities in the ",
  "<em>λ</em> framework are defined for <strong>A</strong><sup><em>k</em></sup>, ",
  "while cohort quantities in the <em>R</em><sub>0</sub> framework are defined for ",
  "<strong>A</strong><sub><em>k</em></sub><sup>&dagger;</sup>. Notice the structural symmetry ",
  "between the <em>λ</em> and <em>R</em><sub>0</sub> frameworks."
)
) |>
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "condensed")
  ) |>
  column_spec(3, width = "7cm") |>
  row_spec(0, extra_css = "border-bottom: 2px solid black;") |>
  footnote(
    general_title = "",
    general = paste0(
      "<div style='border-top: 1px solid black; padding-top: 6px;'>",
      "<span style='font-size: 90%;'>",
      "<strong>Note.</strong> Aggregation with <code>leslie_aggregate()</code> ",
      "changes the projection interval from Δt to <em>k</em>Δt.<br>",
      "<sup>&dagger;</sup> <strong>A</strong><sub><em>k</em></sub><sup></sup> denotes the <em>k</em>-step-ahead projection matrix ",
      "in the <em>R</em><sub>0</sub> framework.</span></div>"
    ),
    footnote_as_chunk = TRUE,
    escape = FALSE
  )

## -----------------------------------------------------------------------------
# New example used for Leslie-to-Leslie aggregation
# These matrices replace the earlier example matrices (A, U, F)
# Leslie population projection matrix for Pimephales promelas (MatrixID = 248187)

# 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 == 248187, ]
# A <- matA(mpm)[[1]]
# U <- matU(mpm)[[1]]
# F <- matF(mpm)[[1]]



# Leslie population projection matrix
A <- rbind(
  c(0.00, 0.00,  6.00, 62.50),
  c(0.50, 0.00,  0.00,  0.00),
  c(0.00, 0.18,  0.00,  0.00),
  c(0.00, 0.00,  0.12,  0.00)
)

# Reproduction matrix
F <- rbind(
  c(0.00, 0.00,  6.00, 62.50),
  c(0.00, 0.00,  0.00,  0.00),
  c(0.00, 0.00,  0.00,  0.00),
  c(0.00, 0.00,  0.00,  0.00)
)

# Survival/transition matrix
U <- rbind(
  c(0.00, 0.00,  0.00,  0.00),
  c(0.50, 0.00,  0.00,  0.00),
  c(0.00, 0.18,  0.00,  0.00),
  c(0.00, 0.00,  0.12,  0.00)
)

# Partitioning matrix used for the aggregation examples below
# (first two age classes combined, last two combined)
# In practice this matrix is constructed automatically inside leslie_aggregate()
# It is reproduced here only to verify the aggregation results.
P <- rbind(
  c(1,1,0,0),
  c(0,0,1,1)
)

# Dimensions used in the aggregation examples
n <- 4     # dimensionality of the original matrix A
m <- 2     # dimensionality of the aggregated matrix
k <- n / m #(also constructed automatically)

## -----------------------------------------------------------------------------

# Leslie-to-Leslie MPM aggregation
# in the lambda framework, standard aggregation
res_leslie_lambda_std <- leslie_aggregate(
  matA = A,
  ngroups = 2,
  framework = "lambda",
  criterion = "standard"
)


# Aggregated Leslie matrix
res_leslie_lambda_std$matAk_agg

# Effectiveness of aggregation
res_leslie_lambda_std$effectiveness

# ---- Check: lambda consistency ----

# Expected lambda
spectral_radius(A)^k

# Lambda of aggregated model
spectral_radius(res_leslie_lambda_std$matAk_agg)

# ---- Check: stable age distribution consistency ----

# Stable age distribution
w <- leslie_stable_age(A=A)

# Expected 
c(P%*%w)

# Aggregated stable age distribution
w_agg <-leslie_stable_age(res_leslie_lambda_std$matAk_agg)
w_agg

# ---- Check: interstage flow consistency ----

# Expected
P%*%(A%^%k)%*%diag(w)%*%t(P)

# Interstage flows from aggregated model
res_leslie_lambda_std$matAk_agg%*%diag(w_agg)

## -----------------------------------------------------------------------------

# Leslie-to-Leslie MPM aggregation
# in the lambda framework, elasticity-consistent aggregation
res_leslie_lambda_el <- leslie_aggregate(
  matA = A,
  ngroups = 2,
  framework = "lambda",
  criterion = "elasticity"
)

# Aggregated Leslie matrix
res_leslie_lambda_el$matAk_agg

# Effectiveness of aggregation
res_leslie_lambda_el$effectiveness

# get aggregated F and U
F_agg <-matrix(0,2,2)
F_agg[1,] <- res_leslie_lambda_el$matAk_agg[1,]
U_agg <- res_leslie_lambda_el$matAk_agg
U_agg[1,] <-0

# ---- Check: lambda consistency ----

# Expected
spectral_radius(A)^k

# Lambda for aggregated model
spectral_radius(res_leslie_lambda_el$matAk_agg)

# ---- Check: stable age distribution consistency ----

# Stable age distribution
w <- leslie_stable_age(A=A)

# Expected
c(P%*%w)

# Stable age distribution of aggregated odel
w_agg <-leslie_stable_age(res_leslie_lambda_el$matAk_agg)
w_agg

# ---- Check: reproductive values consistency ----

# Get reproductive values
v <- leslie_reproductive_values(A=A)

# Expected
c(solve(P%*%diag(w)%*%t(P))%*%P%*%(v*w))

# Reproductive values of aggregated model
leslie_reproductive_values(A=res_leslie_lambda_el$matAk_agg)

# ---- Check: elasticities consistency ----

# Get elasticities of lambda^k
EA <- mpm_elasticity(A%^%k,framework="lambda")$elasticity

# Expected
P%*%EA%*%t(P)

# Elasticities of lambda^k for aggregated model
mpm_elasticity(res_leslie_lambda_el$matAk_agg,framework="lambda")$elasticity

# ---- Check: generation time consistency ----

# Survival/transition matrix
U_k <- U%^%k

# Reproduction matrix
F_k <- A%^%k - U_k

# Expected generation time
generation_time(matF=F_k,matU=U_k,framework="lambda")$generation_time

# Generation time of aggregated model
generation_time(matF=F_agg,matU=U_agg,framework="lambda")$generation_time

## -----------------------------------------------------------------------------

# Leslie-to-Leslie MPM aggregation
# in the R0 framework, standard aggregation
res_leslie_R0_std <- leslie_aggregate(
  matA = A,
  ngroups = 2,
  framework = "R0",
  criterion = "standard"
)

# Aggregated Leslie matrix
res_leslie_R0_std$matAk_agg

# Effectiveness of aggregation
res_leslie_R0_std$effectiveness

# get F and U for aggregated matrix
F_agg <- matrix(0,2,2)
F_agg[1,] <- res_leslie_R0_std$matAk_agg[1,]
U_agg <- res_leslie_R0_std$matAk_agg
U_agg[1,] <- 0

# ---- Check: generation time consistency ----

# Expected
R0 <- spectral_radius(F%*%solve(diag(4)-U))
R0 

# R0 for aggregated model
spectral_radius(F_agg%*%solve(diag(2)-U_agg))


# ---- Check: cohort stable age distribution consistency ----

# Standardized cohort projection matrix
A_ref <- (F + U*R0)/R0

# Standardized cohort projection matrix for aggregated model
A_ref_agg <- (F_agg + U_agg*R0)/R0

# w is the stable age distribution
w <- leslie_stable_age(A=A_ref)

# Expected
c(P%*%w)

# Stable age distribution for aggregated model
w_agg <-leslie_stable_age(A_ref_agg)
w_agg

# ---- Check: cohort interstage flows consistency ----

# Survival/transition matrix (k steps ahead)
U_k <- U%^%k

# Reproduction matrix (k steps ahead)
F_k <- A%^%k - U_k

# A projection matrix (k steps ahead)
A_k <- F_k + U_k

# Expected
P%*%A_k%*%diag(w)%*%t(P)

# Interstage flows for aggregated model
(F_agg + U_agg)%*%diag(w_agg)

## -----------------------------------------------------------------------------

# Leslie-to-Leslie MPM aggregation
# in the R0 framework, elasticity-consistent aggregation
res_leslie_R0_el <- leslie_aggregate(
  matA = A,
  ngroups = 2,
  framework = "R0",
  criterion = "elasticity"
)

# Aggregated Leslie matrix
res_leslie_R0_el$matAk_agg

# Effectiveness of aggregation
res_leslie_R0_el$effectiveness

# Get F and U for aggregated matrix
F_agg <- matrix(0,2,2)
F_agg[1,] <- res_leslie_R0_el$matAk_agg[1,]
U_agg <- res_leslie_R0_el$matAk_agg
U_agg[1,] <- 0

# ---- Check: R0 consistency ----

# Expected
R0 <- spectral_radius(F%*%solve(diag(4)-U))
R0 

# R0 of aggregated model
spectral_radius(F_agg%*%solve(diag(2)-U_agg))


# ---- Check: cohort stage age distribution consistency ----

# Standardized cohort projection matrix
A_ref <- (F + U*R0)/R0

# Standardized cohort projection matrix for aggregated model
A_ref_agg <- (F_agg + U_agg*R0)/R0

# Cohort stable age distribution
w <- leslie_stable_age(A=A_ref)

# Expected
c(P%*%w)

# Cohort stage age distribution for aggregated model
w_agg <-leslie_stable_age(A_ref_agg)
w_agg

# ---- Check: cohort reproductive values consistency ----

# Cohort reproductive values
v <- leslie_reproductive_values(A_ref)

# Expected
c(solve(P%*%diag(w)%*%t(P))%*%P%*%(v*w))

# Cohort reproductive values of aggregated model
leslie_reproductive_values(A_ref_agg)

# ---- Check: elasticities of R0 consistency ----
# These elasticities are normalized to sum to 1

# Survival/transition matrix (k steps ahead)
U_k <- U%^%k

# Reproduction matrix (k steps ahead)
F_k <- A%^%k - U_k

# Cohort projection matrix (k steps ahead)
A_k <- F_k + U_k

# Elasticity of R0 (normalized to sum to 1)
EA <- mpm_elasticity(matF=F_k,matU=U_k,framework="R0")$elasticity

# Expected
P%*%EA%*%t(P)

# Elasticity of R0 for the aggregated model
mpm_elasticity(matF=F_agg,matU=U_agg,framework="R0")$elasticity

# ---- Check: cohort generation time consistency ----

# Expected
generation_time(matF=F_k, matU = U_k, framework="R0")$generation_time

# Generation time of aggregated model
generation_time(matF=F_agg,matU=U_agg, framework="R0")$generation_time

