In petrelharp/landsim: Simulate Populations on Landscapes

fig.dim <- 4
knitr::opts_chunk$set(fig.width=2*fig.dim,fig.height=fig.dim,fig.align='center')
library(Matrix)
library(raster)
library(landsim)

f <- vital( function (N,...) { N+x }, x=3 )
f(10)
f$x
f$x <- 12
f(10)

To run a simulation, we need to know where the species lives, how it disperses, and how it reproduces and dies. Given a demographic model of this sort, we then keep track of the numbers of various genotypes across the species range.

These are divided among different objects as follows:

`population`

Class population:

habitat

: a Raster with values that can be used in other computations.

accessible

: vector of indices of cells in habitat that migrants will attempt to go to

habitable

: vector of indices of cells in habitat that may have positive population

genotypes

: character vector of the genotypes

N

: matrix indexed by (habitable cells) x (genotypes) giving the number of each genotype in each habitable cell

Set-up:

This can be initiated with make_population() by specifying the following:

habitat

: The path to the file where the raster is stored.

inaccessible.value

: The values in the raster that should be marked as inaccessible.

uninhabitable.value

: The values in the raster that should be marked as not habitable.

genotypes

: A character vector of genotypes.

Example:

pop <- make_population(
                 habitat = system.file("extdata/test_raster.gri",package="landsim"),
                 inaccessible.value = "NA",
                 uninhabitable.value = 0.0,
                 genotypes = c( "aa", "aA", "AA" )
             )

`migration`

One structure we use several times is the setup for smoothing (with migrate); for this we need to know:

Class migration:

kern

: function that gives weights for neighboring cells as a function of Euclidean distance

sigma

: scaling factor on distance before being passed to kern

radius

: maximum range of the smoother

normalize

: normalization factor applied to the smoother: the rows of the resulting matrix will sum to this value; if this is NULL no normalization is done.

n.weights

: weights on the number of times to apply the smoother to produce one "migration" step, soif $M$ is the one-step smoother, and the weights are $n=(n_1,n_2,\ldots)$, the resulting smoother is $$ x \mapsto (1-\sum_k n_k) x + \sum_k n_k M^k x . $$

Example:

This is a migration object for a Gaussian kernel with standard deviation 100 and maximum dispersal distance 1000.

migr <- migration(
                 kern = "gaussian",
                 sigma = 100,
                 radius = 1000,
                 normalize = 1.0
             )

and we can modify it so that a migrating individual takes no steps with probability $1/4$, one step with probability $1/2$, or five steps with probability $1/4$:

migr <- migration(
                 kern = "gaussian",
                 sigma = 100,
                 radius = 1000,
                 normalize = 1.0,
                 n.weights = c(1/2,0,0,0,1/4)
             )

migration.matrix

A migration.matrix object extends a migration object, in that it has a precomputed migration matrix, and hence is tied to a particular population setup (although this is not included). This has additional entries

Class migration.matrix:

M

: The pre-computed migration matrix, which has rows and columns indexed by accessible locations in the habitat.

habitable.inds

: The indices in columns of M that correspond to habitable locations.

Example:

This is added with the function setup_migration().

migr <- setup_migration( migr, population=pop )

The migration matrix has rows and columns indexed by accessible locations in the habitat, as opposed to population$N, which is indexed by habitable locations.

`vital_rate`

Vital rates may be numbers or vectors but may also be functions, to allow them to depend on the current state of the system. We would like to have the parameters in these visible, changeable, and carried around with them. To this end, we define vital objects by specifying a function and auxialliary parameters that will be kept along with the function.

Simple example:

This class allows assigning and accessing things directly from the environment of the function, without having them clutter up the global environment.

f <- vital( function (N,...) { N+x }, x=3 )
f(10)
f$x
f$x <- 12
f(10)

Example:

germ.vital <- vital(
             function (N, ...) {
                out <- r0 / ( 1 + migrate(rowSums(N),competition)/carrying.capacity )
                return( cbind( aa=out, aA=s*out, AA=s^2*out ) )
             },
             r0 = 0.01,  # one in ten seeds will germinate at low densities
             s = 1.5,    # multiplicative selective benefit of the A allele
             carrying.capacity = 10,
             competition = migration(
                                     kern="gaussian",
                                     sigma=100,
                                     radius=300,
                                     normalize=1
                                 )
         )

Note that this also has a migration object as part of its definition. This needs to be set up also:

germ.vital <- setup_vital(germ.vital,pop)

`demography`

Several of the base demographic parameters we might want to specify (as a fixed number or a layer), or compute (as a function of the current population state). Such parameters can be either something numeric or else a function that is applied to N to get the result.

To do a single generation we need to know:

Class demography:

prob.seed

: parameter object for probability of seeding per individual per year

fecundity

: parameter object for mean number of seeds per seeding individual per year

prob.germination

: parameter object for probability of germination

prob.survival

: parameter object for probability of survival of already existing individuals

pollen.migration

: migration object for pollen dispersal (normalized to total pollen production)

seed.migration

: migration object for seed dispersal (normalized to 1)

genotypes

: vector of names of the genotypes represented by the layers

mating

: genotype x genotype x genotype tensor, with entry [i,j,k] the probability that genotypes i and j combine to make offspring genotype k

Example:

demog <- demography(
        prob.seed = 0.2,
        fecundity = 100,
        prob.germination = germ.vital,
        prob.survival = 0.9, 
        seed.migration = migr,
        pollen.migration = migration(
                         kern="gaussian",
                         sigma=300,
                         radius=900,
                         normalize=10
                     ),
        genotypes = c("aa", "aA", "AA")
     )