Metapopulation projections
David Garcia-Callejas, Maria Paniw, and cxr team
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_1*N_{s1}
+ \sum_{i \in {s2,...,sn}} \beta_i*Ni + \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.
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"]]
## alpha beta1 beta2 beta3 beta4 beta5
## Sj 1.0388308 1.3457376 0.069977335 -0.078802049 -0.0599899328 0.002809132
## Sn 3.0471897 1.5928148 -0.050277967 -0.044215025 -0.0214356158 -0.012578936
## Sr 3.8510474 2.3649396 -0.037753281 -0.081006020 -0.0408717604 -0.001894939
## Rn 1.2175032 0.0000000 -0.089419233 0.000000000 0.0000000000 0.000000000
## Rr 0.9076192 0.1785048 -0.022440715 -0.029090383 0.0039708918 0.000000000
## D -0.3153730 0.1924553 0.007495063 -0.007089005 0.0001276958 0.000000000
## Ds 0.7000000 0.1200000 0.000000000 0.000000000 0.0000000000 0.000000000
## O 2.6156340 0.7026241 -0.081252277 -0.061627015 -0.0826407224 0.001873481
## beta6 beta7
## Sj 0.004493186 0.0006890698
## Sn -0.020976595 0.0187594251
## Sr -0.052637814 -0.0333495434
## Rn 0.000000000 0.0000000000
## Rr 0.000000000 0.0000000000
## D 0.000000000 0.0000000000
## Ds 0.000000000 0.0000000000
## O -0.026285279 0.0000000000
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]]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0
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:
## [,1] [,2] [,3]
## [1,] 10 10 15
## [2,] 10 10 13
Where there are 10 juvenile individuals at site 1, 15 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"))
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"]]
## alpha beta1 beta2 beta3 beta4 beta5 beta6 beta7
## Sj NA NA NA NA NA NA NA NA
## Sn NA NA NA NA NA NA NA NA
## Sr NA NA NA NA NA NA NA NA
## Rn NA NA NA NA NA NA NA NA
## Rr NA NA NA NA NA NA NA NA
## D NA NA NA NA NA NA NA NA
## Ds NA NA NA NA NA NA NA NA
## O NA NA NA NA NA NA NA NA
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
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0388308163 0.380762792 2.7282887 6.366387e-03
## dens 0.0699773352 0.010531993 6.6442633 3.047375e-11
## env 1.3457375997 0.353523654 3.8066409 1.408671e-04
## densS2 -0.0788020486 0.014226983 -5.5389150 3.043514e-08
## densS3 -0.0599899328 0.027802360 -2.1577281 3.094898e-02
## dens:env 0.0028091316 0.009296323 0.3021766 7.625174e-01
## env:densS2 0.0044931860 0.012676525 0.3544494 7.230022e-01
## env:densS3 0.0006890698 0.023931772 0.0287931 9.770296e-01
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"]]
## alpha beta1 beta2 beta3 beta4 beta5 beta6
## Sj 1.038831 1.345738 0.06997734 -0.07880205 -0.05998993 0.002809132 0.004493186
## Sn NA NA NA NA NA NA NA
## Sr NA NA NA NA NA NA NA
## Rn NA NA NA NA NA NA NA
## Rr NA NA NA NA NA NA NA
## D NA NA NA NA NA NA NA
## Ds NA NA NA NA NA NA NA
## O NA NA NA NA NA NA NA
## beta7
## Sj 0.0006890698
## Sn NA
## Sr NA
## Rn NA
## Rr NA
## D NA
## Ds NA
## O NA
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:
- 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.
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.
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.
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)
