Metapopulation projections
In cxr: A Toolbox for Modelling Species Coexistence in R

knitr::opts_chunk$set(message = FALSE, warning = FALSE)

# library(tidyverse,quietly = TRUE)
# library(Matrix) # needed for fill_demography_matrix

Introduction

The dynamics and composition of natural communities are ultimately determined by how individuals interact with their biotic and abiotic environment via demographic rates (survival and reproduction) in order to persist in space and time. Here, we develop hierarchical functions to scale from stage-specific demographic rates to community dynamics. The functions first describe stage-specific demographic rates for various species that are affected simultaneously by environmental factors and biotic interactions across different habitats. These rates then determine transitions among different life-cycle stages in a local populations and dispersal among populations, which, in turn, characterizes the metapopulation dynamics of all species in a community.

The metapopulation framework implemented in these functions can model the dynamics of any number of species present in any number of sites, with the only constraint that we assume three distinct life stages, which generally correspond to juvenile individuals, non-reproductive adults, and reproductive adults. The methodology follows the vec-permutation approach introduced in Hunter & Caswell (2005) for projecting metapopulation abundances.

The vital rates considered in this framework include the survival probability of each life stage, the transitions between stages, reproductive outputs, and dispersal. Specifically:

Sj: Survival probability of juveniles
Sn: Survival probability of non-reproductive adults
Sr: Survival probability of reproductive adults
Rn: Probability of a non-reproductive adult to transition to the reproductive stage
Rr: Probability of a reproductive adult to remain in the reproductive stage
D: Probability of dispersal (emigration)
Ds: Probability of surviving all stages of dispersal and establishing in a new site (immigration)
O: Recrutiment, i.e., per capita number of offspring that survive to the juvenile stage produced (≥ 0)

These vital rates can be parameterized to differ weakly or strongly between sites and species. They arise from the combined effect of species-specific values (given by the density of the focal species, $N_{s1}$ below), the effect of individuals of other species present in the same site ($N_{s2,...,sn}$), and potentially the effect of environmental variables ($E$). In particular, for a given vital rate $v$ of a species $s1$, they are computed as:

$$Logit(v_{s1}) = \alpha + \beta_1N_{s1} + \sum_{i \in {s2,...,sn}} \beta_iNi + \beta_E*E$$

In this formulation, interactions between species densities and the environment are not shown for readability, but these can be included as well.

# All defaults
knitr::include_graphics("Figure1.png")

Data structures

The vec-permutation approach relies on a block-diagonal formulation of the demographic and dispersal processes. For each species, demographic changes are described by B, a block-diagonal matrix where the $i^{th}$ diagonal block $K_i$ is a 3x3 matrix population model at site $i$. Dispersal is described by $M$, another block-diagonal matrix where the $j^{th}$ diagonal block $L_j$ is a $i*i$ matrix of disperal probabilities for each stage $j$. In our case, non-reproductive adults disperse to become reproductive adults in a different site, and we therefore populate the corresponding off-diagonal element of $M$. Our metapopulation model assumes that local demography happens first and then individuals disperse.

To run our approach of metapopulation dynamics of interacting species, we need to first define a series of constant values and data structures. Note that we explicitly consider environmental forcing: currently, our framework accepts a single environmental factor that varies throughout the number of projected timesteps.

# load the package
library(cxr)

# define species names
sp <- c("s1","s2","s3")
# number of species
num.sp <- length(sp)

# define site names
sites <- c("sa","sb")
# number of sites
num.sites <- length(sites)

# number of demographic stages - this should be always fixed
num.stages <- 3

# vital rates to account for - these names should be fixed
rates <- c("Sj","Sn","Sr","Rn","Rr","D","Ds","O")

# years (or time steps) of simulations
years <- 100 

# simulate environment from normal distribution
set.seed(123)
env <- rnorm(years, mean=0, sd=1)

At this stage, we need to generate or import the coefficients for obtaining vital rates for every species in every site. This data structure thus holds the $\alpha$ and $\beta$ coefficients from Eq. 1. In this example, we import these coefficients in an already created list named metapopulation_example_param. We suggest to separate in different scripts the generation of the data structure and the projection of the dynamics. Nevertheless, at the end of this vignette we explain in more detail the different options to generate or import these coefficients.

data("metapopulation_example_param", package = "cxr")

# coefficients for species "s1" in site "sa"
metapopulation_example_param[["s1"]][["sa"]]

Next, we define the structure that will hold the vec-permutation matrices, with the given number of species, sites, and stages. There are three types of matrices in the vec-permutation methodology: permutation matrices, dispersal, and demography. There is one matrix of each type for each species, and these matrices all have dimensions of number of stages multiplied by number of sites in both rows and columns.

# build the templates for the vec-permutation matrices
# returns a nested list: [[matrix.type]][[species number]]
# where matrix.type is "demography", "dispersal", or "permutation"
vpm <- vec_permutation_matrices(num.sp,num.sites,num.stages)
# for example
vpm[["demography"]][[1]]

We also need to define manually the initial abundances of each species at each stage in each site. These are stored in a list in which each element holds the abundance of a given species in a matrix of dimensions number of sites in rows and life stages in columns.

# initial.densities: list: [[species]][sites*stages]
# this needs to be filled manually
initial.densities <- list()
# sp1
initial.densities[[1]] <- matrix(c(10,10,15,10,10,13),
                                 nrow = num.sites,
                                 ncol = num.stages,
                                 byrow = TRUE)
# sp2
initial.densities[[2]] <- matrix(c(15,5,3,15,5,3),
                                 nrow = num.sites,
                                 ncol = num.stages,
                                 byrow = TRUE)
# sp2
initial.densities[[3]] <- matrix(c(5,5,2,3,4,2),
                                 nrow = num.sites,
                                 ncol = num.stages,
                                 byrow = TRUE)

For example, the initial abundances of the first species are as follows:

initial.densities[[1]]

Where there are r initial.densities[[1]][1,1] juvenile individuals at site 1, r initial.densities[[1]][1,3] reproductive adults at site 1, and so on.

Lastly, we need to define the data structure that will store the whole projected dynamics, i.e. species densities in each timestep, site, and stage. This is, again, a nested list with species in the first level, year (or more generally, timestep) in the second level, and each element of the third level being a matrix of dimensions of number of sites (rows) by number of stages (columns).

# projected.densities list: [[species]][[years]][[sites*stages]
projected.densities <- list()
for(i.sp in 1:num.sp){
  projected.densities[[i.sp]] <- list()
  for(i.year in 1:years){
    projected.densities[[i.sp]][[i.year]] <- matrix(0,
                                                    nrow = num.sites,
                                                    ncol = num.stages)
  }
}

Projecting population densities

The densities of each species at each timestep are calculated within a loop. In this loop, first we update the transition matrix of each species, given their prior densities and coefficients. This is done with the function fill_transition_matrix. Second, we update the vec-permutation demography and dispersal matrices of each species, with the functions fill_demography_matrix and fill_dispersal_matrix. With these two matrices, we can finally obtain the expected density of each species at the next timestep using the function calculate_densities.

# create an auxiliary list, to keep track of the densities per timestep
current.densities <- initial.densities

for(i.year in 1:years){

  # Update projected densities at this timestep
  for(i.sp in 1:num.sp){
    projected.densities[[i.sp]][[i.year]] <- current.densities[[i.sp]]
  }

  # update transition matrices ----------------------------------------------

  # this is a list per species and site
  transition_matrices <- list()

  for(i.sp in 1:length(sp)){
    transition_matrices[[i.sp]] <- list()
    for(i.site in 1:length(sites)){

      # store the transition matrix for this sp and site
      transition_matrices[[i.sp]][[i.site]] <- fill_transition_matrix(focal.sp = i.sp,
                                                                      site = i.site,
                                                                      param = metapopulation_example_param,
                                                                      env = env[i.year],
                                                                      current.densities = current.densities)

    }# for each site
    names(transition_matrices[[i.sp]]) <- sites
  }# for each species
  names(transition_matrices) <- sp

  # update demography and dispersal matrices --------------------------------

  for(i.sp in 1:length(sp)){

    # demography
    vpm[["demography"]][[i.sp]] <- fill_demography_matrix(focal.sp = i.sp,
                                                          vpm = vpm,
                                                          transition_matrices = transition_matrices)
    # dispersal
    vpm[["dispersal"]][[i.sp]] <- fill_dispersal_matrix(focal.sp = i.sp,
                                                        num.sites = num.sites,
                                                        param = metapopulation_example_param,
                                                        vpm = vpm,
                                                        env = env[i.year],
                                                        current.densities = current.densities)
  }# for i.sp

  # update densities --------------------------------------------------------
  # all stages and sites for each species
  for(i.sp in 1:length(sp)){
    current.densities[[i.sp]] <- calculate_densities(focal.sp = i.sp,
                                                     vpm,
                                                     current.densities)
  }# for i.sp
}# for i.year

Finally, we can plot the resulting dynamics over time, species, and/or sites

# transform list of projected densities to dataframe
df <- densities_to_df(projected.densities)

# tidy
df$species <- dplyr::recode(as.character(df$species), "1" = "S1", "2" = "S2", "3" = "S3")

# plot
dynamics.plot <- ggplot2::ggplot(df,ggplot2::aes(year,density,col=species))+
  ggplot2::geom_line()+
  ggplot2::facet_grid(stage~site,scales = "free")+
  ggplot2::scale_color_manual(name="",values=c("darkgreen","orange","darkred"))+
  ggplot2::xlab("Simulation year")+ggplot2::ylab("Total density")+
  ggplot2::theme_bw(base_size=20)+
  ggplot2::theme(panel.grid = ggplot2::element_blank())+
  ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
        panel.grid.minor = ggplot2::element_blank(),
        strip.background = ggplot2::element_blank(),
        panel.border = ggplot2::element_rect(colour = "black"))

dynamics.plot

Appendix: Different options for building the coefficient data structure

Users can generate the list with parameter coefficients through a set of helper functions. The data structure is created with the function build_param, and populated with the function generate_vital_rate_coefs.

The function build_param accepts the name of the species modelled, the number or name of sites, a string giving the vital rates included, and an argument env specifying whether environmental forcing is to be accounted for. Lastly, the coefficients from Eq. 1 may include the interaction between the density of one or more species and the environment. The user has the possibility to limit the number of env:species interactions by setting the argument num.params. For example, given three species and including environmental forcing, there are a minimum of 1+1+3 coefficients (the intercept, one coefficient for the environmental effect, and one coefficient for the effect of each species), and a maximum of 1+1+3+3 coefficients, by including each interaction. In this example, we include all environment:density interactions, for a total of eight coefficients. It is assumed that the interactions included follow the order in which the species are entered.

# this builds an empty data structure
example_param <- build_param(sp = sp,
                     sites = sites,
                     rates = rates,
                     env = env,
                     num.params = 8)

The example_param structure is a nested list, sorted by species in the first level and site in the second level. Each element is a matrix of dimensions number of vital rates by number of coefficients.

example_param[["s1"]][["sa"]]

There are a number of ways in which these matrices can be populated. The elements, as shown by Eq. 1, are equivalent to the coefficients (on the link scale) of a Generalized Linear Model (GLM). Therefore, if the user has a tidy table with the same number of coefficients, this table can be input to the generate_vital_rate_coefs function via the glm.object argument. For example, using an example coefficient table from a GLM:

data("glm_example_coefs",package = "cxr")

glm_example_coefs

Notice that the coefficient names are different from those in the example_param tables. In order to import the GLM table into the structure, we also need to specify the name equivalences:

glm.coef.equivalence <- list("alpha" = "(Intercept)",
                             "beta1" = "env",
                             "beta2" = "dens", 
                             "beta3" = "densS2", 
                             "beta4" = "densS3",
                             "beta5" = "dens:env",
                             "beta6" = "env:densS2",
                             "beta7" = "env:densS3")

And now, we can import these values. Assuming that these coefficients refer to the survival rate of juveniles for species s1 at site sa:

example_param <- generate_vital_rate_coefs(example_param,
                                           sp = "s1",
                                           sites = "sa",
                                           vital.rate = "Sj",
                                           glm.object = glm_example_coefs,
                                           glm.coef.equivalence = glm.coef.equivalence)

# the resulting table. Note that we only entered the coefficients for survival of juveniles (Sj)
example_param[["s1"]][["sa"]]

Otherwise, coefficients can be entered manually:

example_param$s2$sa["Sr",5] <- 0.001

Metapopulaiton example

In metapopulation_example_param, the coefficients for each demographic rate were obtained from GLMs fitted on simulated data from a relatively simple community with interaction across trophic levels: one prey species (S1), one meso-predator (S2), and one top predator (S3). You can think of a system such as rabbit, fox, and lynx.

The underlying data were simulated using the following rules:

knitr::include_graphics("rabbit.png")

(a) The prey (S1) is relatively more sensitive to environmental fluctuations than either the predators. Survival of prey (across all stages) is negatively affected by the predators, and more so by the top predator (S3), which is assumed to be more specialized on the prey.

knitr::include_graphics("fox.png")

(b) The meso-predator (S2) is a generalist, and its survival across all stages is positively affected by prey density, but not as strongly as the survival of the top predator; it competes with the top predator and is more strongly negatively affected by its density than the other way around.

knitr::include_graphics("lynx.png")

(c) The top predator (S3) is specialized on the prey and therefore responds strongly (and positively) to prey abundance. It is negatively affected by the meso-predator, but with caveats as explained above.

(d) All three species disperse between 2 sites. Dispersal is highest for the top predator, followed by the meso-predator, and prey. Dispersal mortality takes the opposite pattern, being highest for the prey.

Generating coefficient randomly

And, finally, coefficients can also be generated by sampling from a normal distribution with a given mean and standard deviation. Below we enter the coefficient for the intercept (vr.coef = "alpha") of the survival of reproductive adults (vital.rate = c("Sr")), for the first species (sp = "s1") and all sites (by not specifying the sites argument), setting the mean and sd values from which random samples will be drawn.

example_param <- generate_vital_rate_coefs(param = example_param,
                                   sp = "s1",
                                   # sites = "sa",
                                   vital.rate = c("Sj"),#,"Sn","Sr","Rn","Rr","D","Ds","O"),
                                   vr.coef = "alpha",
                                   mean.coef = .1,sd.coef = 0)