R/structured_population_growth.R

In gauravsk/ecoevoapps: Simulate Dynamics of Ecology/Evolution Models

#' Simulate growth of a (st)age-structured population
#'
#' @param leslie_mat A leslie matrix
#' @param init A vector of initial abundances in each (st)age
#' @param time Length of time over which to run the simulation
#' @return Returns a matrix with nrow = nrow(leslie_mat) and ncol = time
#' @seealso [plot_structured_population_size()],
#'   [plot_structured_population_agedist()], and
#'   [plot_structured_population_lambda()] for plotting different aspects of the
#'   population trajectory, and [plot_leslie_diagram()] for plotting an
#'   transition diagram based on the specified Leslie matrix
#' @examples
#' leslie_matrix <- matrix(c(0, 8,1, 1, 0.4,0,0,0,0,0.8,0,0,0,0,0.1,0),
#' ncol = 4, byrow = TRUE)
#' structured_pop_init <- c(10,10,10,10)
#' structured_pop_time <- 5
#' run_structured_population_simulation(leslie_mat = leslie_matrix, init =
#' structured_pop_init, time = structured_pop_time)
#' @export
run_structured_population_simulation <- function(leslie_mat = matrix(c(0, 8,1, 1,
                                                                       0.4,0,0,0,
                                                                       0,0.8,0,0,
                                                                       0,0,0.1,0),
                                                                     ncol = 4, byrow = TRUE),
                                           init = c(10,10,10,10),
                                           time = 100) {
  # Check that leslie matrix is a square
  if(!(ncol(leslie_mat) == nrow(leslie_mat))) stop("leslie_mat should be a square matrix")


  pop <- matrix(init, nrow = length(init), ncol = 1) # start the matrix with initial population sizes
  nn <- sum(init) # total population size
  lambda <- numeric() # vector to hold discrete population growth

  # simulate dynamics in a for-loop
  for(i in 1:time) {
    # run one time-step of the model
    pop <- cbind(pop, floor(leslie_mat %*% pop[,i,drop=F]))
    nn[i+1] <- sum(pop[1:3,i+1]) # total population size
    lambda[i] <- nn[i+1]/nn[i] # calculate lambda for the whole population
    if(nn[i+1]==0) { break } # stop if population crashes
    # stop when stages reached equilibrium
    if(i > 10) {
      if(round(lambda[i],3) == round(lambda[i-2],3)){
        break
      }
    }
  }
  rownames(pop) <- paste0("Age ", 1:nrow(pop))
  return(pop)
}

#' Converts population matrix into a long-format dataframe
#'
#' For internal use only
#' @param pop_growth_matrix matrix where each row is an age class and each column
#' is a time step
#' @import dplyr
#' @keywords internal
popmat_to_df <- function(pop_growth_matrix) {
  # To suppress CMD Check
  time <- NULL

  popgrowth_df <-
    pop_growth_matrix %>%
    t() %>%
    as_tibble %>%
    dplyr::mutate(time = row_number()) %>%
    tidyr::pivot_longer(cols = -time, names_to = "Age class",
                 values_to = "Population size")
  return(popgrowth_df)
}

#' Generate a trajectory of population size in each (st)age over time
#'
#' @param pop_growth_matrix population growth matrix, as generated by
#'   `run_structured_population_simulation` (each row is an age class and each
#'   column is a time step)
#' @return ggplot object of population trajectory over time
#' @seealso [run_structured_population_simulation()] to simulate the growth of a
#'   structured population given a Leslie matrix,
#'   [plot_structured_population_agedist()] and
#'   [plot_structured_population_lambda()] for plotting different aspects of the
#'   population trajectory, and [plot_leslie_diagram()] for plotting an
#'   transition diagram based on the specified Leslie matrix
#' @examples
#' leslie_matrix <- matrix(c(0, 8,1, 1, 0.4,0,0,0,0,0.8,0,0,0,0,0.1,0),
#' ncol = 4, byrow = TRUE)
#' structured_pop_init <- c(10,10,10,10)
#' structured_pop_time <- 50
#' structured_pop_out <- run_structured_population_simulation(leslie_mat = leslie_matrix, init =
#' structured_pop_init, time = structured_pop_time)
#' plot_structured_population_size(structured_pop_out)
#' @import ggplot2
#' @export
plot_structured_population_size <- function(pop_growth_matrix) {
  # To suppress CMD Check
  time <- `Population size` <- `Age class` <- NULL

  pop_growth_df <- popmat_to_df(pop_growth_matrix)
  # if(eigen(leslie)$values[1]<1) {tmp = .8} else {tmp = .2}

  trajs <-
    ggplot(pop_growth_df, aes(x = time, y = `Population size`, group = `Age class`)) +
    geom_line(aes(linetype = `Age class`, col = `Age class`)) +
    labs(title = "Population growth simulation") +
    xlab("Time step") +
    ylab("Number of individuals") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_log10(expand = c(0, 0)) +
    ecoevoapps::theme_apps()

  return(trajs)
}

#' Generate a trajectory of the proportion of individuals in each (st)age over
#' time
#'
#' @param pop_growth_matrix population growth matrix, as generated by
#'   `run_structured_population_simulation` (each row is an age class and each
#'   column is a time step)
#' @param leslie_mat Leslie matrix used to generate the population trajectory
#'   (optional). If the Leslie matrix is provided, the trajectory also includes
#'   lines indicating the stable age distribution
#' @seealso [run_structured_population_simulation()] to simulate the growth of a
#'   structured population given a Leslie matrix,
#'   [plot_structured_population_size()] and
#'   [plot_structured_population_lambda()] for plotting different aspects of the
#'   population trajectory, and [plot_leslie_diagram()] for plotting an
#'   transition diagram based on the specified Leslie matrix
#' @examples
#' leslie_matrix <- matrix(c(0, 8,1, 1, 0.4,0,0,0,0,0.8,0,0,0,0,0.1,0),
#' ncol = 4, byrow = TRUE)
#' structured_pop_init <- c(10,10,10,10)
#' structured_pop_time <- 50
#' structured_pop_out <- run_structured_population_simulation(leslie_mat = leslie_matrix, init =
#' structured_pop_init, time = structured_pop_time)
#' plot_structured_population_agedist(structured_pop_out, leslie_matrix)
#' @return ggplot of age distribution over time
#' @import ggplot2
#' @export
plot_structured_population_agedist <- function(pop_growth_matrix, leslie_mat = NULL) {
  # To suppress CMD Check
  `Population size` <- time <- prop <- `Age class` <- NULL


  pop_growth_df_prop <- popmat_to_df(pop_growth_matrix) %>%
    group_by(factor(time)) %>%
    dplyr::mutate(prop = `Population size`/sum(`Population size`))
  # if(eigen(leslie)$values[1]<1) {tmp = .8} else {tmp = .2}

  trajs <-
    ggplot(pop_growth_df_prop, aes(x = time, y = prop, group = `Age class`)) +
    geom_line(aes(linetype = `Age class`, col = `Age class`)) +
    labs(title = "Trajectory of age structure ") +
    xlab("Time step") +
    ylab("Proportion of individuals") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    theme_apps()

  # Add stable age distribution from Leslie matrix
  # if the leslie matrix is provided to the function
  if(!is.null(leslie_mat)) {
    trajs <-
      trajs +
      geom_hline(yintercept=c(Re(eigen(leslie_mat)$vec[,1])/
                                sum(Re(eigen(leslie_mat)$vec[,1]))),
                 linetype="dashed", show.legend=F) +
      labs(caption = "Dashed black lines show stable age distribution
           calculated from the first eigenvector of the Leslie matrix") +
      theme(plot.caption = element_text(size = 8))

  }
  return(trajs)
}


#' Generate a trajectory of the total population size over time
#'
#' @param pop_growth_matrix population growth matrix, as generated by
#' `run_structured_population_simulation` (each row is an age class and each column
#' is a time step)
#' @param leslie_mat Leslie matrix used to generate the population trajectory (optional).
#' If the Leslie matrix is provided, the trajectory also includes lines indicating
#' the stable age distribution
#' @import ggplot2
#' @importFrom stats na.omit
#' @seealso [run_structured_population_simulation()] to simulate the growth of a
#'   structured population given a Leslie matrix,
#'   [plot_structured_population_size()] and
#'   [plot_structured_population_agedist()] for plotting different aspects of the
#'   population trajectory, and [plot_leslie_diagram()] for plotting an transition
#' @examples
#' leslie_matrix <- matrix(c(0, 8,1, 1, 0.4,0,0,0,0,0.8,0,0,0,0,0.1,0),
#' ncol = 4, byrow = TRUE)
#' structured_pop_init <- c(10,10,10,10)
#' structured_pop_time <- 50
#' structured_pop_out <- run_structured_population_simulation(leslie_mat = leslie_matrix, init =
#' structured_pop_init, time = structured_pop_time)
#' plot_structured_population_lambda(structured_pop_out, leslie_matrix)
#' @export
plot_structured_population_lambda <- function(pop_growth_matrix,
                                              leslie_mat = NULL) {
  # To suppress CMD Check
  total_popsize <- time <- lambda <- NULL


  pop_growth_df <-
    t(pop_growth_matrix) %>%
    as_tibble() %>%
    dplyr::mutate(total_popsize = rowSums(.),
           lambda = total_popsize/lag(total_popsize),
           time = row_number()) %>%
    na.omit() # eliminate the first row, as there is no lambda for time step 1

  lambda_plot <-
    ggplot(pop_growth_df, aes(x = time, y = lambda)) +
    geom_line(color = "grey50") +
    labs(title = "Overall population growth rate") +
    xlab("Time step") +
    ylab(expression("Population growth rate (" ~lambda* ")")) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    theme_apps()


  # Add stable age distribution from Leslie matrix
  # if the leslie matrix is provided to the function
  if(!is.null(leslie_mat)) {
    lambda_plot <-
      lambda_plot +
      geom_hline(yintercept = Re(eigen(leslie_mat)$val[1])) +
      labs(caption = "Solid black lines shows the long-term population growth
           rate calculated from the dominant eigenvalue of the Leslie matrix") +
      theme(plot.caption = element_text(size = 8))

  }
  return(lambda_plot)
}

#' Generate diagram of population structure
#'
#' @param leslie_mat Leslie matrix from which to generate diagram
#' @import diagram
#' @seealso [run_structured_population_simulation()] to simulate the growth of a
#'   structured population given a Leslie matrix,
#'   [plot_structured_population_size()],
#'   [plot_structured_population_agedist()], and
#'   [plot_structured_population_lambda()] for plotting different aspects of the
#'   population trajectory
#' @examples
#' leslie_matrix <- matrix(c(0, 8,1, 1, 0.4,0,0,0,0,0.8,0,0,0,0,0.1,0),
#' ncol = 4, byrow = TRUE)
#' plot_leslie_diagram(leslie_matrix)
#' @export
plot_leslie_diagram <- function(leslie_mat) {
  n_ages <- nrow(leslie_mat)
  name_vec <- parse(text = paste0("Age[",1:n_ages,"]"))
  diagr <- plotmat(A = leslie_mat, pos = n_ages,
                            curve = 0.5, lwd = 1.5, my = -0.1,
               name = name_vec,
               arr.length = 0.2, arr.width = 0.25, arr.lwd = 2,
               arr.type = "simple", self.lwd = 2, self.shiftx = 0.115,
               self.shifty = 0.1, self.cex = .5, box.size = 0.1,
               dtext = 0.2, box.lwd = 3.5, main="Life cycle diagram"
               )

  return(diagr)

}