Skip to contents

Mechanistic Metapopulation Simulator

Source: R/sim.R
sim.Rd

This function simulates population growth and dispersal providing a given environmental scenario. All parameters for the simulation must be set in advance using initialise.

Usage

sim(
  obj,
  time,
  burn = 0,
  return_mu = FALSE,
  cl = NULL,
  progress_bar = FALSE,
  quiet = TRUE
)

Arguments

obj

sim_data object created by initialise containing all simulation parameters and necessary data

time

positive integer vector of length 1; number of time steps simulated

burn

positive integer vector of length 1; the number of burn-in time steps that are discarded from the output

return_mu

logical vector of length 1; if TRUE demographic process return expected values; if FALSE the rpois function should be used

cl

an optional cluster object created by makeCluster needed for parallel calculations

progress_bar

logical vector of length 1 determines if progress bar for simulation should be displayed

quiet

logical vector of length 1; determines if warnings should be displayed

Value

This function returns an object of class sim_results which is a list containing the following components:

  • extinction - TRUE if population is extinct or FALSE otherwise

  • simulated_time - number of simulated time steps without the burn-in ones

  • N_map - 3-dimensional array representing spatio-temporal variability in population numbers. The first two dimensions correspond to the spatial aspect of the output and the third dimension represents time.

In case of a total extinction, a simulation is stopped before reaching the specified no. of time steps. If the population died out before reaching the burn threshold, then nothing can be returned and an error occurs.

Details

This is the main simulation module. It takes the sim_data object prepared by initialise and runs simulation for a given number of time steps. The initial (specified by the burn parameter) time steps are skipped and discarded from the output. Computations can be done in parallel if the name of a cluster created by makeCluster is provided.

Generally, at each time step, simulation consists of two phases: local dynamics and dispersal.

Local dynamics (which connects habitat patches in time) is defined by the function growth. This parameter is specified while creating the obj using initialise, but can be later modified by using the update function. Population growth can be either exponential or density-dependent, and the regulation is implemented by the use of Gompertz or Ricker models (with a possibility of providing any other, user defined function). For every cell, the expected population density during the next time step is calculated from the corresponding number of individuals currently present in this cell, and the actual number of individuals is set by drawing a random number from the Poisson distribution using this expected value. This procedure introduces a realistic randomness, however additional levels of random variability can be incorporated by providing parameters of both demographic and environmental stochasticity while specifying the sim_data object using the initialise function (parameters r_sd and K_sd, respectively).

To simulate dispersal (which connects habitat patches in space), for each individual in a given cell, a dispersal distance is randomly drawn from the dispersal kernel density function. Then, each individual is translocated to a randomly chosen cell at this distance apart from the current location. For more details, see the disp function.

Examples

if (FALSE) {

# data preparation
library(terra)

n1_small <- rast(system.file("input_maps/n1_small.tif", package = "rangr"))
K_small <- rast(system.file("input_maps/K_small.tif", package = "rangr"))

sim_data <- initialise(
  n1_map = n1_small,
  K_map = K_small,
  r = log(2),
  rate = 1 / 1e3
)

# simulation
sim_1 <- sim(obj = sim_data, time = 100)

# simulation with burned time steps and progress bar
sim_2 <- sim(obj = sim_data, time = 100, burn = 20, progress_bar = TRUE)

# example with parallelization
library(parallel)
cl <- makeCluster(detectCores())

# parallelized simulation
sim_3 <- sim(obj = sim_data, time = 100, cl = cl, progress_bar = TRUE)
stopCluster(cl)


# visualisation
plot(n1_small)

plot(to_rast(
  sim_1,
  time_points = 100,
  template = sim_data$K_map
))

plot(to_rast(
  sim_1,
  time_points = c(1, 10, 50, 100),
  template = sim_data$K_map
), range = range(sim_1$N_map, na.rm = TRUE))

}