Simulate population dynamics through time.

Use functions from anadrofish to simulate population change through time relative to upstream and downstream passage probabilities and uncertainty in life-history information.

Usage

sim_pop(
  species = c("ALE", "AMS", "BBH"),
  nyears = 50,
  river,
  max_age = NULL,
  nM = NULL,
  fM = 0,
  n_init = runif(1, 1e+06, 8e+08),
  spawnRecruit = NULL,
  eggs = NULL,
  sr = 0.5,
  b = 0.21904,
  s_juvenile = NULL,
  upstream = 1,
  downstream = 1,
  downstream_j = 1,
  output_years = c("last", "all"),
  age_structured_output = FALSE,
  sex_specific = TRUE,
  custom_habitat = NULL
)

Arguments

species: Species for which population dynamics will be simulated. Choices include American shad ("AMS"), alewife ("ALE"), and blueback herring ("BBH").
nyears: Number of years for simulation.
river: River basin. Available rivers implemented in package can be viewed by calling get_rivers with no arguments (e.g., get_rivers()). Alternatively, the user can specify rivers = sample(get_rivers, 1) to randomly sample river within larger simulation studies. Information about each river can be found in the habitat dataset.
max_age: Maximum age of fish in population. If NULL (default), then based on the maximum age of females for the corresponding region in the max_ages dataset.
nM: Instantaneous natural mortality. If NULL (default), then based on the average of males and females for the corresponding region in the mortality dataset.
fM: Instantaneous fishing mortality. The default value is zero.
n_init: Initial population seed (number of Age-1 individuals) used to simulate the starting population. Default is to use a draw from a wide uniform distribution, but it may be beneficial to narrow once expectations for abundance at population stability are determined.
spawnRecruit: Probability of recruitment to spawn at age. If NULL (default), then probabilities are based on the mean of male and female recruitment to first spawn at age from the maturity dataset.
eggs: Number of eggs per female. Can be a vector of length 1 if eggs per female is age invariant, or can be vector of length max_age if age-specific. If NULL (default) then estimated based on weight-batch fecundity regression relationships for each life-history region in the olney_mcbride dataset (Olney and McBride 2003) and mean number of batches spawned (6.1 +/- 2.1, McBride et al. 2016) using make_eggs for American shad or parameters from Sullivan et al. (2019) for river herring.
sr: Sex ratio (expressed as percent female or P(female)).
b: Density-dependent parameter for the Beverton-Holt stock-recruitment relationship. The default value (0.21904) is used to approximate a larval carrying capacity at 100 adult fish per acre based on fecundity American shad and river herring. A value of about 0.05 approximates a carrying capacity that corresponds to about 500 adult fish per acre.
s_juvenile: Survival from hatch to outmigrant. If NULL (default) then simulated from a (log) normal distribution using mean and sd of Sc through 70 d for 1979-1982 from Crecco et al. (1983) in crecco_1983.
upstream: Numeric of length 1 representing proportional upstream passage through dams.
downstream: Numeric of length 1 indicating proportional downstream survival through dams.
downstream_j: Numeric of length 1 indicating proportional downstream survival through dams for juveniles.
output_years: Temporal level of detail provided in output. The default value of 'last' returns the final year of simulation. Any value other than the default 'last' will return data for all years of simulation. This is useful for testing.
age_structured_output: Should population and spawner abundance in the output dataframe be age-structured? If FALSE (default), then pop (non-spawning population) and spawners (spawning population) are summed across all age classes for each year of simulation. If TRUE then pop and spawners are returned for each age class. For the sake of managing outputs, abundances for pop and spawners are reported for all age classes 1-13 regardless of max_age, but all abundances for ages greater than max_age are zero.
sex_specific: Whether to use sex-specific life-history data.
custom_habitat: A dataframe containing columns corresponding to the those in the output from custom_habitat_template. The default, NULL uses the default habitat data set for a given combination of species and river.

Value

A dataframe containing simulation inputs (arguments to sim_pop) and outputs (number of spawners) by year, minimally including the following variables:

river Name of river.
region Regional grouping, see get_region.
govt Governmental unit at downstream terminus of habitat unit.
lat Latitude at downstream terminus of habitat unit.
habitat Amount of accessible habitat given upstream passage and river.
year year of simulation. If output_years = "last" then nyears.
upstream upstream passage through dams. Currently a single value.
downstream adult downstream survival through dams. Currently a single value.
downstream_j adult downstream survival through dams. Currently a single value.
max_age_m maximum age of males in population.
max_age_m maximum age of females in population.
nM_m instantaneous natural mortality of males.
nM_f instantaneous natural mortality of females.
fM instantaneous fishing mortality rate in population.
n_init number of adult fish used to seed population.
sr sex ratio of adults.
s_juvenile survival of juveniles from hatch to outmigrant.
s_spawn_m survival of males during upstream spawning migration.
s_spawn_f survival of females during upstream spawning migration.
s_postspawn_m survival of males after spawning.
s_postspawn_f survival of females after spawning.
iteroparity probability of repeat spawning.
spawners number of spawning adults entering river each year.
pop total population abundance prior to spawning each year.
juveniles_out number of juveniles exiting river each year.

Additional columns and years are included depending on the values passed to ouput_years, age_structured_output.

References

Sullivan, K.M, M.M. Bailey, and D.L. Berlinksky. 2019. Digital Image Analysis as a Technique for Alewife Fecundity Estimation in a New Hampshire River. North American Journal of Fisheries Management 39:353-361.

McBride, R. S., R. Ferreri, E. K. Towle, J. M. Boucher, and G. Basilone, G. 2016. Yolked oocyte dynamics support agreement between determinate- and indeterminate-method estimates of annual fecundity for a northeastern United States population of American shad. PLoS ONE 11(10):10.1371/journal.pone.0164203

Olney, J. E. and R. S. McBride. 2003. Intraspecific variation in batch fecundity of American shad (*Alosa sapidissima*): revisiting the paradigm of reciprocal latitudinal trends in reproductive traits. American Fisheries Society Symposium 35:185-192.

Examples

# Example usage
if (FALSE) { # \dontrun{
  
# Example 1. Simulate a single population one time for 50 years ----------------
  library(anadrofish)
  library(ggplot2)
  
  # Set rng seed ----
  set.seed(1)
  
  # Simulate a single population for fifty years one time
  res <- sim_pop(
    species = "AMS",
    nyears = 50,
    river = 'Susquehanna',
    fM = 0,
    n_init = runif(1, 10e5, 50e6),
    sr = rbeta(1, 100, 100),
    upstream = 1,
    downstream = 1,
    downstream_j = 1,
    output_years = 'all',
    age_structured_output = FALSE)
  
  # Plot the output - the result will be different for each
  # simulated model run above
  ggplot(res, aes(x = year, y = spawners)) +
    geom_line() +
    xlab("Year") +
    ylab("Number of spawners")
   

# Example 2. Simulate a single population 100 times in parallel ----------------
  # Package load ----
  library(snowfall)
  library(parallel)
  library(anadrofish)
  library(tidyverse)
  
  # Set rng seed ----
  set.seed(1)
  
  # Parallel settings ----
  # Get number of cores
  ncpus = parallel::detectCores() - 1
  
  # Initialize snowfall
  sfInit(parallel = TRUE, cpus = ncpus, type="SOCK")
  
  # Wrapper function ----
  sim <- function(x){

  # . Call simulation ----
  # Run with a single set of upstream and downstream
  # dam passage probabilities
  res <- sim_pop(
    species = "AMS",
    nyears = 50,
    river = 'Susquehanna',
    fM = 0,
    n_init = runif(1, 8e6, 50e6),
    sr = rbeta(1, 100, 100),
    upstream = 1,
    downstream = 1,
    downstream_j = 1,
    output_years = 'all',
    age_structured_output = FALSE)
  
  # . Define the output lists ----
  retlist <- list(res=res)      
  
  return(retlist)    
  
  }  
    
  # Parallel execution ----
  # . Load libraries on workers -----
  sfLibrary(anadrofish)
  
  # . Distribute to workers -----
  # Number of simulations to run
  niterations <- 1000
  
  # . Run the simulation ----
  start <- Sys.time()
  
  result <- sfLapply(1:niterations, sim) 
  
  total_time <- Sys.time() - start
  total_time
  
  # . Stop snowfall ----
  sfStop()
  
  # Results ----
  # 'result' is a list of lists. Save this:
  # save(result, file = "sim_result.rda")
  
  # Extract results dataframes by string and rbind them
  res <- lapply(result, function(x) x[[c('res')]])
  library(data.table)
  resdf <- data.frame(rbindlist(res))
  
  glimpse(resdf)
  
  # Results summary
  plotter <- resdf %>% 
    group_by(year) %>% 
    summarize(fit = mean(spawners),
              lwr = quantile(spawners, 0.025),
              upr = quantile(spawners, 0.975))
    
  # Plotting
  ggplot(plotter, aes(x = year, y = fit)) +
    geom_line() +
    geom_ribbon(aes(xmax = year, ymin = lwr, ymax = upr,
                    color = NULL), alpha = 0.40)
    
  
# Example 3. Multi-river simulation --------------------------------------------
# Simulate population size for randomly
# selected rivers and randomly chosen passage
# probabilities from a pre-defined list.
    
# Run simulations in parallel, saving age-structured
# output, but only for the final year of simulation
  # Package load ----
  library(snowfall)
  library(parallel)
  library(anadrofish)
  library(tidyverse)
  
  # Set rng seed ----
  set.seed(1)
  
  # Parallel settings ----
  # Get number of cores
  ncpus = parallel::detectCores() - 1
  
  # Initialize snowfall
  sfInit(parallel = TRUE, cpus = ncpus, type="SOCK")
  
  # Wrapper function ----
  sim <- function(x){
  
    # . Call simulation ----
    # Define passage scenarios (ASFMC 2020) 
    passages <- list(
      c(0,0,0),
      c(1,1,1),
      c(0.4, 0.80, 0.95))
    passage <- unlist(sample(passages, 1))  
      
    # Run the sim  
    res <- sim_pop(
      species = "AMS",
      nyears = 50,
      river = get_rivers(species = "AMS")[sample(1:length(get_rivers(species = "AMS")), 1)],    
      fM = 0,
      n_init = runif(1, 10e5, 50e6),
      sr = rbeta(1, 100, 100),
      upstream = passage[1],
      downstream = passage[2],
      downstream_j = passage[3],
      output_years = 'last',
      age_structured_output = TRUE
      )
    
    # . Define the output lists ----
    retlist <- list(res = res)      
  
    return(retlist)      
  }  
  
  # Parallel execution ----
  # . Load libraries on workers -----
  sfLibrary(anadrofish)
  
  # . Distribute to workers -----
  # Number of simulations to run
  niterations <- 1000
  
  # . Run the simulation ----
  start <- Sys.time()
  
  result <- sfLapply(1:niterations, sim) 
  
  total_time <- Sys.time() - start
  total_time
  
  # . Stop snowfall ----
  sfStop()
  
  # Results ----
  # 'result' is a list of lists. Save this:
  # save(result, file = "sim_result.rda")
  
  # Extract results dataframes by string and rbind them
  res <- lapply(result, function(x) x[[c('res')]])
  library(data.table)
  resdf <- rbindlist(res)
  
  # . Summary statistics by passage scenario -----
  resdf$river <- as.character(resdf$river)
  
  # Sum population size and spawners across age groups
  resdf$spawners <- resdf %>%
    select(grep("spawners", colnames(resdf))) %>%
    rowSums 
  resdf$pop <- resdf %>%
    select(grep("pop", colnames(resdf))) %>%
    rowSums 
  
  # Assign scenarios based on upstream passage
  resdf$scenario <- 'No dams'
  resdf$scenario[resdf$upstream==0] <- 'No passage'
  resdf$scenario[resdf$upstream==.40] <- 'Current condition'
  
  # . System-specific summaries ----
  rivers <- resdf %>% 
    group_by(region, river, scenario) %>% 
    summarize(fit = mean(spawners),
              lwr = quantile(spawners, 0.025),
              upr = quantile(spawners, 0.975))
  

  # . Coastal summary ----
  coastal <- rivers %>% 
    group_by(scenario) %>% 
    summarize(fit = sum(fit),
              lwr = sum(lwr),
              upr = sum(upr))
  
  # . River plot ----
  rivers %>% 
    filter(region == "SI") %>% 
    ggplot(aes(x = river, y = fit, color = scenario, fill = scenario)) +
      geom_col() +
      coord_flip() +
      xlab("River") +
      ylab("Number of spawning adults") +
      ggtitle("Southern Iteroparous Region") +
      theme(plot.title = element_text(hjust = 0.5))
  
  # . Coastal plot ----
  ggplot(coastal, aes(x = scenario, y = fit)) +
    geom_point(size = 4) +
    geom_linerange(aes(xmax = scenario, ymin = lwr, ymax = upr),
                   linewidth = 1) +
    xlab("Scenario") +
    ylab("Number of spawning adults")

} # }