A Simple Example With Migration
In CKMRpop: Forward-in-Time Simulation and Tallying of Pairwise Relationships

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

In this vignette we show how to use CKMRpop/spip in a multi-population version with asymmetrical migration rates. It would be good to read the "about_spip" vignette beforehand to understand how migration is implemented in spip.

Our goal here is to model three "populations" of a fish species in which migration occurs during the larval stage, and it is influenced by currents that flow down the coast and make migration occur primarily in one direction. Each of the populations will have the same life tables, etc. We will make it simple.

First load some packages:

library(tidyverse)
library(CKMRpop)

Basic Demography

The basic demography for each, is given below. We assume no mortality from newborn to age one. In effect we are only modeling the the newborns that survive to age 1. Only 3, 4, and 5 year-olds mate, and older fish have somewhat more reproductive success.

pars <- list()
pars$`max-age` <- 5
pars$`fem-surv-probs` <- c(1, 0.7, 0.8, 0.8, 0.8)
pars$`male-surv-probs` <- c(1, 0.7, 0.8, 0.8, 0.8)
pars$`fem-prob-repro` <- c(0, 0, 1, 1, 1)
pars$`male-prob-repro` <- c(0, 0, 1, 1, 1)
pars$`fem-asrf` <- c(0, 0, .5, .7, 1)
pars$`male-asrp` <- c(0, 0, .5, .7, 1)
pars$`fem-rep-disp-par` <- 0.25
pars$`male-rep-disp-par` <- 0.25
pars$`sex-ratio` <- 0.5

We will run these simulations for 25 years and just assume an annual cohort size of 200 fish in each population, and we will go for a constant population size and start off with the stable age distribution.

pars$`number-of-years` <- 25

# given cohort sizes of 250 the stable age distribution
# can be found and used for the initial number of indivs
L <- leslie_from_spip(pars, 250)


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

pars$`cohort-size` <- "const 250"

Sampling

We will simulate a sampling program in which juveniles are sampled when they recruit after their larval phase. We assume that the larval phase starts at birth, but is over by age 1. So we can sample new recruits via prekill sampling at age 1. We will pretend the sampling is non-lethal. We will also non-lethally sample a fraction of the older age-class fish.

pars$`discard-all` <- 0
pars$`gtyp-ppn-fem-pre` <- "10-25 0.05 0.00 0.00 0.00 0.00"
pars$`gtyp-ppn-male-pre` <- "10-25 0.05 0.00 0.00 0.00 0.00"
pars$`gtyp-ppn-fem-post` <- "10-25 0.00 0.01 0.01 0.01 0.01"
pars$`gtyp-ppn-male-post` <- "10-25 0.00 0.01 0.01 0.01 0.01"

OK! At this point, we have set up the demography for a single population. Let's just show that we could simulate a single population with that demography and the sampling like this:

set.seed(15)
one_pop_dir <- run_spip(pars = pars, num_pops = 1)

But, of course, the question is, how do we simulate multiple, interacting populations?

Simulating Multiple Populations

To simulate multiple populations in spip, from within CKMRpop, we set the num_pops option of run_spip() to the number of populations. And then, for the pars option we pass a list of num_pops different lists, each one holding the demography and sampling information for each of the populations in the simulation. In our case, since we are simulating three populations with the same demography, num_pops = 3 and we can construct our pars argument like this:

pars_list <- list(
  pars,
  pars,
  pars
)

because each population has the same demography.

Note that we have not included any migration between the the populations yet, but we can still run a simulation with all of them like this:

# NOTE the following if()...else() blocks are here
# to test whether spip has been installed yet.
# If spip is not available (for example, on CRAN's build machines) this
# is noted and stored package data are used for the variable
# "slurped" to build the remainder of the vignette.
if(spip_exists()) {
  message("spip is installed and will be used")
  set.seed(15)
  three_pop_dir <- run_spip(
    pars = pars_list,
    num_pops = 3
  )
  slurp3 <- slurp_spip(three_pop_dir, num_generations = 1)
} else {
  message("Using stored package data because spip is not installed")  
  slurp3 <- three_pops_no_mig_slurped_results
}

set.seed(15)
three_pop_dir <- run_spip(
  pars = pars_list,
  num_pops = 3
)
slurp3 <- slurp_spip(three_pop_dir, num_generations = 2)

And you can investigate that output by making the plots and other summaries shown in the main vignette. Many of these summaries account for the multiple populations like this:

ggplot_census_by_year_age_sex(slurp3$census_postkill)

Those are the postkill census sizes in the three different populations, which get the labels 0, 1, and 2, according to the order in which they were passed to spip (i.e., their order in the pars_list, above.)

Migration between multiple populations

The above simulated multiple populations, but did not simulate any migration between them.

Here, we will imagine different rates of self-recruitment of the larvae, which means different rates at which surviving larvae are going to be migrants. The populations are labeled 0, 1, and 2, and, as stated above, we will simulate somewhat unidirectional migration between them, like so:

20% of pop 0 individuals are migrants; 80% of those end up in pop 1 and 20% in pop 2
15% of pop 1 indivs are migrants; 95% of them end up in pop 2 and 5% in pop 0
5% of pop 2 are migrants; 90% of those end up in pop 1 and 10% make it to pop 0.

The way we implement this is by adding options to each population's demography parameters like so:

# out-migration rates pops 0, 1, 2
pars_list[[1]]$`fem-prob-mig-out` <- "5-25 1 .20"
pars_list[[1]]$`male-prob-mig-out` <- "5-25 1 .20"

pars_list[[2]]$`fem-prob-mig-out` <- "5-25 1 .15"
pars_list[[2]]$`male-prob-mig-out` <- "5-25 1 .15"

pars_list[[3]]$`fem-prob-mig-out` <- "5-25 1 .05"
pars_list[[3]]$`male-prob-mig-out` <- "5-25 1 .05"

# in-migration rates
pars_list[[1]]$`fem-prob-mig-in` <-  "5-25 1 0.00 0.80 0.20"
pars_list[[1]]$`male-prob-mig-in` <- "5-25 1 0.00 0.80 0.20"

pars_list[[2]]$`fem-prob-mig-in` <-  "5-25 1 0.05 0.00 0.95"
pars_list[[2]]$`male-prob-mig-in` <- "5-25 1 0.05 0.00 0.95"

pars_list[[3]]$`fem-prob-mig-in` <-  "5-25 1 0.10 0.90 0.00"
pars_list[[3]]$`male-prob-mig-in` <- "5-25 1 0.10 0.90 0.00"

Now, we can simulate this:

# NOTE the following if()...else() blocks are here
# to test whether spip has been installed yet.
# If spip is not available (for example, on CRAN's build machines) this
# is noted and stored package data are used for the variable
# "slurped" to build the remainder of the vignette.
if(spip_exists()) {
  message("spip is installed and will be used")
  set.seed(15)
  mig_dir <- run_spip(
    pars = pars_list,
    num_pops = 3
  )
  slurp_mig <- slurp_spip(mig_dir, num_generations = 1)
} else {
  message("Using stored package data because spip is not installed")  
  slurp_mig <- three_pops_with_mig_slurped_results
}

set.seed(15)
mig_dir <- run_spip(
  pars = pars_list,
  num_pops = 3
)
slurp_mig <- slurp_spip(mig_dir, num_generations = 1)

Now, the object slurp_mig holds information about census sizes of the populations and the relationships between samples, etc. That can be broken down by populations.

First, a record of every migration event is stored in slurp_mig$migrants. This can be easily summarized:

slurp_mig$migrants %>%
  count(age, from_pop, to_pop)

For a further example, to tabulate parent-offspring pairs found between different populations, we can compile the relationships and summarize:

# compile relationships
crel <- compile_related_pairs(slurp_mig$samples)

# count number of PO pairs by which populations
# the members were born in
crel %>%
  filter(dom_relat == "PO") %>%
  mutate(
    parent_born_pop = case_when(
      upper_member == 1 ~ born_pop_1,
      upper_member == 2 ~ born_pop_2,
      TRUE ~ NA_integer_
    ),
    child_born_pop = case_when(
      upper_member == 1 ~ born_pop_2,
      upper_member == 2 ~ born_pop_1,
      TRUE ~ NA_integer_
    )
  ) %>% 
  count(parent_born_pop, child_born_pop)