Simulation from species 1 life history

library(tidyverse)
library(CKMRpop)

For this first example, we use the hypothetical life history of species 1. First we have to set spip up to run with that life history.

Setting the spip parameters

spip has a large number of demographic parameters. Typically spip is run as a command-line program in Unix. In CKMRpop, all that action goes on under the hood, but you still have to use the spip parameters. This vignette is not about using spip. For a short listing of all the spip options, do this:

library(CKMRpop)
spip_help()

If you want a full, complete, long listing of all the spip options, then you can do:

library(CKMRpop)
spip_help_full()

All of the “long-form” options to spip are given on the Unix command line starting with two dashes, like --fem-surv-probs. To set parameters within CKMRpop to send to spip, you simply make a named list of input values. The names of the items in the list are the long-format option names without the leading two dashes. For an example, see the package data object species_1_life_history, as described below.

Basic life history parameters

These parameters are included in the package in the variable species_1_life_history. It is named list of parameters to send to spip. The list names are the names of the spip options. It looks like this:

species_1_life_history
#> $`max-age`
#> [1] 20
#> 
#> $`fem-surv-probs`
#>  [1] 0.75 0.76 0.76 0.77 0.77 0.78 0.78 0.79 0.79 0.80 0.80 0.80 0.81 0.81 0.82
#> [16] 0.82 0.82 0.82 0.82 0.82
#> 
#> $`male-surv-probs`
#>  [1] 0.75 0.76 0.76 0.77 0.77 0.78 0.78 0.79 0.79 0.80 0.80 0.80 0.81 0.81 0.82
#> [16] 0.82 0.82 0.82 0.82 0.82
#> 
#> $`fem-prob-repro`
#>  [1] 0.00 0.00 0.00 0.02 0.09 0.36 0.75 0.94 0.99 1.00 1.00 1.00 1.00 1.00 1.00
#> [16] 1.00 1.00 1.00 1.00 1.00
#> 
#> $`male-prob-repro`
#>  [1] 0.00 0.00 0.00 0.00 0.01 0.05 0.22 0.64 0.91 0.98 1.00 1.00 1.00 1.00 1.00
#> [16] 1.00 1.00 1.00 1.00 1.00
#> 
#> $`fem-asrf`
#>  [1] 0.00 0.00 0.00 3.56 3.72 3.88 4.04 4.20 4.36 4.52 4.68 4.84 5.00 5.16 5.32
#> [16] 5.48 5.64 5.80 5.96 6.12
#> 
#> $`male-asrp`
#>  [1] 0.00 0.00 0.00 3.56 3.72 3.88 4.04 4.20 4.36 4.52 4.68 4.84 5.00 5.16 5.32
#> [16] 5.48 5.64 5.80 5.96 6.12
#> 
#> $`offsp-dsn`
#> [1] "negbin"
#> 
#> $`fem-rep-disp-par`
#> [1] 0.7
#> 
#> $`male-rep-disp-par`
#> [1] 0.7
#> 
#> $`mate-fidelity`
#> [1] 0.3
#> 
#> $`sex-ratio`
#> [1] 0.5

We want to add instructions to those, telling spip how long to run the simulation, and what the initial census sizes should be.

So, first, we copy species_1_life_history to a new variable, SPD:

SPD <- species_1_life_history

Now, we can add things to SPD.

Setting Initial Census, New Fish per Year, and Length of Simulation

The number of new fish added each year is called the “cohort-size”. Once we know that, we can figure out what the stable age distribution would be given the survival rates, and we can use that as our starting point. There is a function in the package that helps with that:

# before we tell spip what the cohort sizes are, we need to 
# tell it how long we will be running the simulation
SPD$`number-of-years` <- 100  # run the sim forward for 100 years

# this is our cohort size
cohort_size <- 300

# Do some matrix algebra to compute starting values from the
# stable age distribution:
L <- leslie_from_spip(SPD, cohort_size)

# then we add those to the spip parameters
SPD$`initial-males` <- floor(L$stable_age_distro_fem)
SPD$`initial-females` <- floor(L$stable_age_distro_male)

# tell spip to use the cohort size
SPD$`cohort-size` <- paste("const", cohort_size, collapse = " ")

Specifying the fraction of sampled fish, and in different years

Spip let’s you specify what fraction of fish of different ages should be sampled in different years. Here we do something simple, and instruct spip to sample 1% of the fish of ages 1, 2, and 3 (after the episode of death, see the spip vignette…) every year from year 50 to 75.

samp_frac <- 0.03
samp_start_year <- 50
samp_stop_year <- 75
SPD$`discard-all` <- 0
SPD$`gtyp-ppn-fem-post` <- paste(
  samp_start_year, "-", samp_stop_year, " ", 
  samp_frac, " ", samp_frac, " ", samp_frac, " ",
  paste(rep(0, SPD$`max-age` - 3), collapse = " "),
  sep = ""
)
SPD$`gtyp-ppn-male-post` <- SPD$`gtyp-ppn-fem-post`

Running spip and slurping up the results

There are two function that do all this for you. The function run_spip() runs spip in a temporary directory. After running spip, it also processes the output with a few shell scripts. The function returns the path to the temporary directory. You pass that temporary directory path into the function slurp_spip() to read the output back into R. It looks like this:

set.seed(5)  # set a seed for reproducibility of results
spip_dir <- run_spip(pars = SPD)  # run spip
slurped <- slurp_spip(spip_dir, 2) # read the spip output into R

Note that setting the seed allows you to get the same results from spip. If you don’t set the seed, that is fine. spip will be seeded by the next two integers in current random number sequence.

If you are doing multiple runs and you want them to be different, you should make sure that you don’t inadvertently set the seed to be the same each time.

Some functions to summarize the runs

Although during massive production simulations, you might not go back to every run and summarize it to see what it looks like, when you are parameterizing demographic simulations you will want to be able to quickly look at observed demographic rates and things. There are a few functions in CKMRpop that make this quick and easy to do.

Plot the age-specific census sizes over time

This is just a convenience function to make a pretty plot so you can check to see what the population demographics look like:

ggplot_census_by_year_age_sex(slurped$census_postkill)

This shows that the function leslie_from_spip() does a good job of finding the initial population numbers that accord with the stable age distribution.

Assess the observed survival rates

We can compute the survival rates like this:

surv_rates <- summarize_survival_from_census(slurped$census_postkill)

That returns a list. One part of the list is a tibble with observed survival fractions. The first 40 rows look like this:

surv_rates$survival_tibble %>%
  slice(1:40)
#> # A tibble: 40 x 7
#>     year   pop   age sex        n cohort surv_fract
#>    <int> <int> <int> <chr>  <int>  <int>      <dbl>
#>  1    20     0    20 female     1      0        0  
#>  2    20     0    19 female     2      1        0.5
#>  3    21     0    20 female     1      1        0  
#>  4    20     0    18 female     1      2        1  
#>  5    21     0    19 female     1      2        1  
#>  6    22     0    20 female     1      2        0  
#>  7    20     0    17 female     2      3        1  
#>  8    21     0    18 female     2      3        1  
#>  9    22     0    19 female     2      3        1  
#> 10    23     0    20 female     2      3        0  
#> # … with 30 more rows

The second part of the list holds a plot with histograms of age-specific, observed survival rates across all years. The blue line is the mean over all years.

surv_rates$plot_histos_by_age_and_sex

To compare these values to the parameter values for the simulation, you must pass those to the function:

surv_rates2 <- summarize_survival_from_census(
  census = slurped$census_prekill,
  fem_surv_probs = SPD$`fem-surv-probs`, 
  male_surv_probs = SPD$`male-surv-probs`
)

# print the plot
surv_rates2$plot_histos_by_age_and_sex

Here, the red dashed line is the value chosen as the parameter for the simulations. The means are particularly different for the older age classes, which makes sense because there the total number of individuals in each of those year classes is smaller.

The distribution of offspring number

It makes sense to check that your simulation is delivering a reasonable distribution of offspring per year. This is the number of offspring that survive to just before the first prekill census. Keep in mind that, for super high-fecundity species, we won’t model every single larva, we just don’t start “keeping track of them” until they reach a stage that is recognizable in some way.

We make this summary from the pedigree information. In order to get the number of adults that were present, but did not produce any offspring, we also need to pass in the postkill census information. Also, to get lifetime reproductive output, we need to know how old individuals were when they died, so we also pass in the information about deaths.

To make all the summaries, we do:

offs_and_mates <- summarize_offspring_and_mate_numbers(
  census_postkill = slurped$census_postkill,
  pedigree = slurped$pedigree,
  deaths = slurped$deaths, lifetime_hexbin_width = c(1, 2)
)

Note that we are setting the lifetime reproductive output hexbin width to be suitable for this example.

The function above returns a list of plots, as follows:

Age and sex specific number of offspring

offs_and_mates$plot_age_specific_number_of_offspring

Especially when dealing with viviparous species (like sharks and mammals) it is worth checking this to make sure that there aren’t some females having far too many offspring.

Lifetime reproductive output as a function of age at death

Especially with long-lived organisms, it can be instructive to see how lifetime reproductive output varies with age at death.

offs_and_mates$plot_lifetime_output_vs_age_at_death

Yep, many individuals have no offspring, and you have more kids if you live longer.

Fractional contribution of each age class to the year’s offspring

Out of all the offspring born each year, we can tabulate the fraction that were born to males (or females) of each age. This summary shows a histogram of those values. The represent the distribution of the fractional contribution of each age group each year.

offs_and_mates$plot_fraction_of_offspring_from_each_age_class

The blue vertical lines show the means over all years.

The distribution of the number of mates

Some of the parameters in spip affect the distribution of the number of mates that each individual will have. We can have a quick look at whether the distribution of number of mates (that produced at least one offspring) appears to be what we might hope it to be.

mates <- count_and_plot_mate_distribution(slurped$pedigree)

That gives us a tibble with a summary, like this:

head(mates$mate_counts)
#> # A tibble: 6 x 6
#>   sex     year   pop parent   num_offs num_mates
#>   <chr>  <int> <int> <chr>       <int>     <int>
#> 1 female    20     0 F0_0_0          2         1
#> 2 female    20     0 F10_0_1         5         4
#> 3 female    20     0 F10_0_11        3         1
#> 4 female    20     0 F10_0_12        1         1
#> 5 female    20     0 F10_0_13        2         1
#> 6 female    20     0 F10_0_2         3         2

And also a plot:

mates$plot_mate_counts