Population introduction via translocation functions

In poems: Pattern-Oriented Ensemble Modeling System

knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  fig.show = "asis",
  fig.align = "centre",
  results = "markup",
  comment = "#>"
)

In this vignette we present a simple example of the poems workflow to show how to use a custom translocation function to introduce populations at defined time steps. The vignette also shows how to define different growth rates for populations in different regions.

Note that the examples are for a single stage matrix model.

Setup

As before, We begin by loading the necessary packages and setting a temporary output directory.

library(poems)
library(raster)
library(sf)
library(scales)
library(stringi) # for randomly generating file names.

OUTPUT_DIR <- tempdir()

# function to round to any arbitrary value
round_any <- function(x, accuracy, f = round) {
  f(x/ accuracy) * accuracy
}

\newpage

Step 1: Build the population model for the study region

Create a model template using the PopulationModel class. This model is spatially explicit, so is defined using the Region class. In addition we make the model temporally explicit so we can introduce populations at defined time steps.

Study region

First, we'll define our study region. For this example, we utilise a raster::RasterLayer() of Thylacine habitat suitability.

# Region raster
tasmania_raster <- poems::tasmania_raster
tasmania_raster

# Equal area projection
tasPrj <- 'PROJCS["Tasmania_Lambert_Azimuthal",
                 GEOGCS["GCS_WGS_1984",
                        DATUM["D_WGS_1984",
                              SPHEROID["WGS_1984",6378137.0,298.257223563]],
                        PRIMEM["Greenwich",0.0],
                        UNIT["Degree",0.0174532925199433]],
                 PROJECTION["Lambert_Azimuthal_Equal_Area"],
                 PARAMETER["False_Easting",0.0],
                 PARAMETER["False_Northing",0.0],
                 PARAMETER["Central_Meridian",147],
                 PARAMETER["Latitude_Of_Origin",-42.2],
                 UNIT["Meter",1.0]]'

# Template raster to project to
tempExt <- projectExtent(tasmania_raster, tasPrj)
res(tempExt) <- 10000 # 10 km resolution
tempExt

# Project the region
tasmania_raster <- projectRaster(tasmania_raster, tempExt, 
                                 method = "ngb")
plot(tasmania_raster, main = "Tasmania raster", 
     legend = FALSE,
     col = "#2E8B57", colNA = "grey75")

\newpage

Now we can define our poems::Region

# Tasmania study region (735 non-NA cells stored in the order shown) #
region <- Region$new(template_raster = tasmania_raster)
region$region_raster

# Establish HS template and starting location #
# This will be our initial introduction point
int_ll <- sf_project(from = "EPSG:4326",
                     to = tasPrj,
                     pts = cbind(146.44, -41.18))
int_point <- region$region_indices[
  which(region$region_indices == 
          cellFromXY(tasmania_raster, xy = int_ll))]

# row which corresponds to initial introduction site
int_index <- which(region$region_indices == int_point) # 114

# plot of region, and introduction locations
plot(region$region_raster, 
     main = "Tasmanian study region (cell indices)",
     colNA = "grey75",
     addfun = function() {
       points(xyFromCell(region$region_raster, int_point), col = "red", pch = 16)
     })

Land-use modifier

Here we read in a land-use modifier layer which we can use to make our region spatiotemporally explicit. This has the effect of altering the HS values through time which causes dynamic changes in habitat suitability, and thus population abundances.

While we could change these values now and supply them as is to the simulator. In this example, we will use a Generator object later on to apply the HS scaling.

# read in the land-use modifier
hs_mod <- poems::tasmania_modifier

plot(hs_mod, zlim = c(0, 1), colNA = "grey75",
     col = hcl.colors(100, "RdYlGn"))

# Habitat suitability
hs_raster <- poems::thylacine_hs_raster
hs_raster <- projectRaster(hs_raster, region$region_raster, method = "bilinear")
hs_raster <- stretch(hs_raster, minv = 0, maxv = 1)
hs_raster

# initial_hs needed for generator
initial_hs <- hs_raster <- stack(replicate(n = nlayers(hs_mod), hs_raster))

\newpage

While we have read in the land-use modifier, we have not applied it to the HS values yet. This means that our HS values are still static as shown by the plot below.

# static HS for the moment
## capacity generator will make it temporally dynamic
plot(hs_raster, zlim = c(0, 1), colNA = "grey75",
     col = hcl.colors(100, "RdYlGn"),
     addfun = function() {
       points(xyFromCell(region$region_raster, int_point), pch = 16, cex = 0.5)
     })

\newpage

Environmental correlation

Next, we'll define a distance-based spatial correlation for applying environmental stochasticity within our model. The generated correlation data is compacted for computational efficiency (with large-scale models).

# Distance-based environmental correlation (via a compacted Cholesky decomposition)
env_corr <- SpatialCorrelation$new(region = region, 
                                   amplitude = 0.496, 
                                   breadth = 80, 
                                   distance_scale = 1000)
env_corr$calculate_compact_decomposition(decimals = 4)

Spatially-varying growth rates

Here we use IBRA regions to define varying growth rates for different populations. These could be expected to occur for example in a wide-ranging species that occurs in varying habitats (e.g. red foxes in Australia).

# allow growth rates to vary by region using IBRA regions
# Tasmania study Interim Bioregionalisation of Australia (IBRA) bioregion cell distribution
ibra_raster <- poems::tasmania_ibra_raster
ibra_raster <- projectRaster(ibra_raster, region$region_raster, method = "ngb")
plot(ibra_raster, colNA = "grey75",
     breaks = seq(1, 9, 1),
     main = "IBRA regions of Tasmania",
     col = hcl.colors(10, "Lajolla"))

ibra_data <- poems::tasmania_ibra_data
ibra_data

# Calculate cell indices and counts for IBRA bioregions
ibra_indices <- lapply(as.list(ibra_data$index),
                       function(i) {
                         which(ibra_raster[region$region_indices] == i)
                       })
str(ibra_indices)

ibra_polygons <- rasterToPolygons(ibra_raster, dissolve = TRUE, na.rm = TRUE)
ibra_polygons@data <- merge(ibra_polygons@data, ibra_data, 
                            by.x = "layer", by.y = "index")
ibra_polygons

plot(ibra_polygons, col = hcl.colors(9, "Lajolla"), border = "black")
text(ibra_polygons, labels = "abbr", cex = 1.2, halo = TRUE)

rmax_regional <- ibra_raster

# seed is set to keep example results constant
{set.seed(27); rmax <- round(rlnorm(9, 0.94, 0.3), 1)}
for (val in 1:9) {
  rmax_regional[rmax_regional == val] <- rmax[val]
}
plot(rmax_regional, colNA = "grey75",
     legend = FALSE, main = "regional growth rates",
     zlim = range(rmax),
     addfun = function() {
       plot(ibra_polygons, border = hcl.colors(9, "Lajolla"), 
            col = NA, add = TRUE)
       text(ibra_polygons, labels = rmax, halo = TRUE)
     },
     col = hcl.colors(100, "Zissou"))

# set upper and lower growth rates per region
ibra_rmax <- cbind(ibra_data, rmax_lower = round(rmax * 0.6, 2),
                   rmax_mean = round(rmax, 2),
                   rmax_upper = round(rmax / 0.75, 2))
ibra_rmax

\newpage

Step 2: Setup the translocation function

Here we generate a custom Translocation class object. This object can be used to translocate populations from one location to another, or as shown here, to introduce populations at defined locations and timesteps. The function could also be expanded to introduce varying number of animals at each site. For simplicity sake, this example releases the same number of animals at all locations.

Introduction sites and times

Here we define the introduction sites by matching region$region_indices to our introduction locations.

We also define the timesteps that the introductions should occur. These timesteps are sequential from 1, and are not defined by, for example, calendar years.

# set up translocation locations and order
intro_trans_ll <- sf_project(from = "EPSG:4326",
                       to = tasPrj,
                       pts = cbind(c(148.01, 144.7, 147.9, 148.27, 145.24), 
                                   c(-40.8, -40.7, -43.2, -42.02, -42.3)))
intro_trans_ll
intro_trans_point <- region$region_indices[which(region$region_indices %in% 
                                             cellFromXY(region$region_raster, 
                                                        xy = intro_trans_ll))]
intro_trans_point <- intro_trans_point[-1]
intro_cells <- intro_trans_point; intro_cells

intro_times <- c(2, 3, 6, 8)

# Introduction times and locations
cbind(intro_times, intro_cells)

plot(region$region_raster,
     main = "Introduction sites",
     col = hcl.colors(100, "Lajolla"),
     addfun = function() {
       plot(ibra_polygons, border = "black", col = NA, add = TRUE)
       points(xyFromCell(region$region_raster, intro_cells), pch = 16, 
              cex = 1.5, col = c("darkgreen", "blue2", "black", "goldenrod"))
       points(region$coordinates[which(region$region_indices %in% intro_cells), ],
       col = "firebrick", cex = 1.5, lwd = 2)
       })

\newpage

Translocation function

Here we define the custom translocation function. It's simply a list object, with a nested function that uses parameters from the model to make changes to the simulated populations.

# User-defined translocation function (list-nested) and alias ####

translocation <- list(

  # Function parameters (passed to function in params list)
  intro_cells = intro_cells, # cells where pops are introduced
  intro_timesteps = intro_times, # timesteps when introduced
  trans_n = 50, # translocated abundances. If not provided by LHS == 50
  region_indices = region$region_indices,

  # Function definition
  translocation_function = function(params) {
    # Unpack parameters (used at every time step)
    intro_cells <- params$intro_cells
    intro_timesteps <- params$intro_timesteps
    simulator <- params$simulator
    stages <- params$stages
    populations <- params$populations
    abundances <- params$abundance
    region_indices <- params$region_indices
    tm <- params$tm # timestep
    sa <- params$stage_abundance
    trans_n <- params$trans_n
    # if introduction at timestep, introduce pops
    if (tm %in% intro_timesteps) {
      # take stage abundance at timestep
      new_sa <- array(sa, c(stages, populations))
      # identifies location of introduction
      trans_loc <- which(region_indices == intro_cells[which(intro_timesteps == tm)])
      # add n individuals regardless of K
      new_sa[trans_loc] <-  new_sa[trans_loc] + trans_n
      return(new_sa)
    } else {
      # else return pops as they are
      new_sa <- array(sa, c(stages, populations))
      return(new_sa)
    }
  }
)
translocation_aliases <- list(intro_cells = "translocation$intro_cells",
                              intro_times = "translocation$intro_timesteps",
                              trans_n = "translocation$trans_n",
                              region_indices = "translocation$region_indices")

\newpage

Step 3: Build generators for dynamically generating model parameters

Growth rate generator

Now we can build a Generator class object that will generate random growth rates based on quantiles between the min and max in each IBRA region. This allows us to sample across the range of values within each region.

# Build a Rmax generator based on sampled IBRA Rmax range quantile
rmax_gen <- Generator$new(
  description = "Rmax",
  spatial_correlation = env_corr,
  generate_rasters = FALSE,
  ibra_data_rmax = ibra_rmax,
  ibra_indices = ibra_indices,
  region_cells = region$region_cells,
  inputs = c("rmax_quantile"),
  outputs = c("growth_rate_max"),
  generative_requirements = list(growth_rate_max = "function"))

# growth_rate_max template
rmax_gen$add_function_template(
  "growth_rate_max",
  function_def = function (params) {
    growth_rate_max <- array(0, params$region_cells)
    for (i in 1:nrow(params$ibra_data_rmax)) {
      growth_rate_max[params$ibra_indices[[i]]] <- 
        stats::qunif(params$rmax_quantile, 
                     min = params$ibra_data_rmax$rmax_lower[i],
                     max = params$ibra_data_rmax$rmax_upper[i])
    }
    return(growth_rate_max)
  },
  call_params = c("ibra_data_rmax", "ibra_indices", "region_cells", "rmax_quantile"))

# test rmax generator at median values
rmax_gen_ex <- rmax_gen$generate(input_values = list(rmax_quantile = 0.5))
rmax_regional[region$region_indices] <- rmax_gen_ex$growth_rate_max
plot(rmax_regional, main = "median regional rmax",
     col = hcl.colors(100),
     addfun = function() {
       plot(ibra_polygons, border = "black", col = NA, add = TRUE)
     })

The map above shows the result of sampling our possible rmax values to their median values. The plot below shows that we can sample across a range of values using the same generator. This functionality allows us to pass rmax_quantile as a variable in a LatinHypercubeSampler class.

# Test multiple quantiles
test_rmax <- lapply(seq(0, 1, 0.1), function(i) {
  region$raster_from_values(rmax_gen$generate(input_values = list(rmax_quantile = i))$growth_rate_max)
})
test_rmax <- stack(test_rmax)
names(test_rmax) <- paste0("Q", seq(0,1,0.1))

# plot
plot(test_rmax, colNA = "grey75",
     legend = TRUE,
     zlim = c(min(values(test_rmax), na.rm = TRUE),
              max(values(test_rmax), na.rm = TRUE)),
     addfun = function() {
       plot(ibra_polygons, border = "black", col = NA, add = TRUE)
     },
     col = hcl.colors(100))

\newpage

Dispersal generator

This generator controls dispersal across our region based on mean dispersal distance and proportion of dispersers.

# Dispersal generator ####
# Set for veriable mean distance, max hard-coded at 150
dispersal_gen <- DispersalGenerator$new(
  region = region,
  spatial_correlation = env_corr,
  generate_rasters = FALSE,
  dispersal_max_distance = 150,
  distance_classes = seq(5, 150, by = 10),
  distance_scale = 1000, # in km
  dispersal_friction = DispersalFriction$new(), # modify coastline distances
  inputs = c("dispersal_p", "dispersal_b"), # proportion and average distance
  decimals = 4)
dispersal_gen$calculate_distance_data()
head(dispersal_gen$distance_data$base, 10)
table(dispersal_gen$distance_data$base$distance_class)

Like with the rmax_gen defined above, we can test our dispersal generator over a range of values to make sure the dispersal curve looks appropriate. We can see in the plot below that as the proportion of dispersal decreases, so does the maximum dispersal distance.

# plot dispersal curves for mean dispersal rates
disp_fun <- function(p, b, distance) {
  p * exp(-distance/b)
}

disp_mat <- data.frame(p = round(runif(1000, 5, 40)/100,2), # prop
                       b = round(runif(1000, 5, 40)) # mean distance
                       )
head(disp_mat)

disp_test <- lapply(1:nrow(disp_mat), function(i) {
  p <- disp_mat[i, "p"]
  b <- disp_mat[i, "b"]
  disp_x <- disp_fun(p, b, seq(5, 150, 5))
  return(disp_x)
})

{par(mar = c(4,4,0.5,0.5))
  matplot(x = seq(5, 150, 5), y = rep(NA, 30), type = "l", ylim = c(0, 0.4),
        xlab = "Disp. dist (km)", ylab = "Prop. disp.", yaxt = "n", xaxt = "n")
  axis(1, at = seq(0, 150, 10))
  axis(2, at = seq(0, 40, 5)/100, labels = seq(0, 40, 5))
  lapply(disp_test, function(i) {
    matplot(x = seq(5, 150, 5), y = unlist(i), type = "l", add = TRUE,
            col = c("#C9C9C944"))}
    )
  lines(x = seq(5, 150, 5),
        y = apply(as.data.frame(disp_test), 1, mean), col = "firebrick")
  }
dev.off()

Now that we're happy the curve follows the right form, we can generate some dispersal data to test our population model with.

# Generate sampled dispersals for p = 0.35, b = 40 (km)
sample_dispersal_data <- dispersal_gen$generate(
  input_values = list(dispersal_p = 0.35, dispersal_b = 40))$dispersal_data
head(sample_dispersal_data[[1]], 10) # examine

\newpage

Capacity generator

Here we build a carrying capacity generator. Carrying capacity is based on maximum density values and scaled by HS (i.e. cells with a HS of 1 contain the highest densities).

The capacity generator is set-up in such a way that it requires multiple parameters to work.

1.max_dens: the maximum theoretical density of populations

2.q_thresh: the quantile threshold used to rescale the HS values

3.trans_n: the number of animals that are introduced. This value is consistent with the translocation function.

capacity_gen <- Generator$new(
  description = "capacity",
  spatial_correlation = env_corr,
  generate_rasters = FALSE,
  time_steps = ncol(initial_hs),
  hs_raster = initial_hs[region$region_indices], # provide full stack of HS. Template attached
  hs_mod = hs_mod[region$region_indices], # provide full stack of LULC modifier. Template attached
  int_index = int_point,
  trans_n = translocation$trans_n, # number of animals introduced
  region_indices = region$region_indices,
  inputs = c("max_dens", "q_thresh", "trans_n"),
  outputs = c("initial_abundance", "carrying_capacity"),
  generative_requirements = list(initial_abundance = "function",
                                 carrying_capacity = "function"))

capacity_gen$add_function_template(
  param = "initial_abundance",
  function_def = function (params) {
    distr_a <- params$hs_raster[, 1]
    ## 0 everywhere except the intro point at the first time step
    ## intro point to trans_n
    ## Could be above or below carrying capacity
    idx <- which(params$region_indices == params$int_index)
    distr_a[idx] <- params$trans_n
    distr_a[-idx] <- 0
    return(distr_a)
  },
  call_params = c("hs_raster", "int_index", "region_indices", "trans_n"))

capacity_gen$add_function_template(
  "carrying_capacity",
  function_def = function (params) {
    idx <- which(params$region_indices == params$int_index)
    distr_k <- params$hs_raster
    distr_mod <- params$hs_mod
    stopifnot("hs_raster and hs_mod have different number of layers" =
                dim(distr_k) == dim(distr_mod))
    # stretch HS values based on q_thresh
    distr_k <- scales::rescale(distr_k, from = c(0, params$q_thresh), to = c(0, 1))
    distr_k[distr_k < 0] <- 0
    distr_k[distr_k > 1] <- 1
    # multiply thresholded HS by hs_modifier
    distr_k <- distr_k * distr_mod
    # rescale back to {0, 1}
    qMax <- max(distr_k, na.rm = TRUE)
    distr_k <- scales::rescale(distr_k, from = c(0, qMax), to = c(0, 1))
    distr_k[distr_k < 0] <- 0
    distr_k[distr_k > 1] <- 1
    # carrying capacity = (HS * maximum density)
    distr_k <- ceiling(distr_k * params$max_dens)
    distr_k[idx, 1] <- params$max_dens
    # distr_k[-idx, 1] <- 0
    return(distr_k)
  },
  call_params = c("hs_raster", "hs_mod", "int_index", "region_indices", "max_dens", "q_thresh"))

# have all parameters been specified correctly
capacity_gen$generative_requirements_satisfied()

Now we have defined our generator we can make some test data for the model

# Generate example initial abundance and declining carrying capacity time-series
generated_k <- capacity_gen$generate(input_values = list(max_dens = 100, q_thresh = 0.90, 
                                                         trans_n = 60))
example_initial_abundance <- generated_k$initial_abundance
example_carrying_capacity <- generated_k$carrying_capacity

# Plot the example initial abundance
example_initial_n_raster <- region$raster_from_values(example_initial_abundance)
example_initial_n_raster
plot(example_initial_n_raster, main = "Example initial abundance",
     col = hcl.colors(100, "Lajolla", rev = TRUE), colNA = "grey75",
     addfun = function() {
       plot(ibra_polygons, border = "black", col = NA, add = TRUE)
     })

# Plot the carrying capacity
## carrying capacity is forced to maximum theoretical value at first time step
example_k <- region$raster_from_values(example_carrying_capacity)
example_k[[c(1, 6, 11)]]
plot(example_k,
     col = hcl.colors(100, "RdYlGn", rev = TRUE), colNA = "grey75",
     addfun = function() {
       plot(ibra_polygons, border = "black", col = NA, add = TRUE)
     },
     zlim = c(0, 100))

\newpage

Step 4: Build a template model

Using the generators we've built we can now test if our simple population model works as expected.

# Template model ####
model_template <- PopulationModel$new(
  region = region,
  time_steps = 11,
  years_per_step = 1,
  stage_matrix = 1, # single-stage
  populations = region$region_cells, # 735
  demographic_stochasticity = TRUE,
  standard_deviation = 0.18,
  density_dependence = "logistic", # Ricker
  harvest = FALSE, # No harvest
  dispersal = dispersal_gen,
  translocation = translocation,
  dispersal_source_n_k = list(threshold = 0.92, cutoff = 0),
  simulation_order = c("translocation", "results", "transition", "dispersal"),
  random_seed = 20230210,
  attribute_aliases = translocation_aliases,
  results_selection = c("abundance"))

model <- model_template$clone()
model$set_attributes(initial_abundance = example_initial_abundance,
                     carrying_capacity = example_carrying_capacity,
                     growth_rate_max =  rmax_gen_ex$growth_rate_max,
                     translocation = translocation,
                     trans_n = 75, # passed through to translocation function
                     dispersal = sample_dispersal_data)
# run poems simulator
results <- population_simulator(model)
results$all$abundance

# timeseries of total abundance
plot(1:11, results$all$abundance, type = "l", 
     xlab = "timestep", ylab = "Total abundance")

The template model runs successfully and now we can make some maps of the dispersal patterns of the populations from their initial introduction.

abund_ras <- region$raster_from_values(results$abundance)
abund_ras[[c(1, 6, 11)]]
abd_max <- round_any(max(values(abund_ras), na.rm = TRUE), 20, f = ceiling)

# plot of abundances. log(x+1) transformed.
plot(log1p(abund_ras),
     col = hcl.colors(100),
     colNA = "grey75",
     addfun = function() {
       plot(ibra_polygons, border = "black", col = NA, add = TRUE)
     },
     zlim = c(0, log1p(abd_max)))

To make sure our translocation function is working correctly, we can run the template model again but with the translocation function turned off.

model$set_attributes(initial_abundance = example_initial_abundance,
                     carrying_capacity = example_carrying_capacity,
                     growth_rate_max = rmax_gen_ex$growth_rate_max,
                     translocation = NULL,
                     dispersal = sample_dispersal_data)
results_notransn <- population_simulator(model) # run poems simulator
results_notransn$all$abundance
results$all$abundance
abund_ras_notransn <- region$raster_from_values(results_notransn$abundance)
abund_ras_notransn[[c(1, 6, 11)]]
diff_ras <- abund_ras - abund_ras_notransn
diff_ras[[1:2]]

We can see from the diff_ras above that the difference between timesteps 2 of the two model runs is 75 individuals, which was the value passed through to the translocation function of the original model run.

And now we can plot the difference.

plotmax <- round_any(max(abs(values(diff_ras)), na.rm = TRUE), 10, ceiling)
plot(diff_ras[[c(1, 6, 11)]],
     zlim = c(-plotmax, plotmax),
     breaks = c(-plotmax, -100, -50, -20, 0, 20, 50, 100, plotmax),
     col = hcl.colors(9, "PuOr"),
     colNA = "grey75",
     addfun = function() {
       plot(ibra_polygons, col = NA, add = TRUE)
     })

\newpage

Step 5: Run multiple simulations

In order to explore the model parameter space to find the best models, we generate Latin hypercube samples of model and generator parameters to be simulated, using the LatinHypercubeSampler class. This class has functionality for generating sample parameters via Uniform, Normal, Lognormal, Beta, and Triangular distributions. For our example we only generate 10 samples. Typically however, a user would need to generate thousands to tens of thousands, of samples.

Define the latin-hypercube for sampling

# Latin-hypercube sampler ####
lhs_gen <- LatinHypercubeSampler$new()

# Habitat suitability threshold
lhs_gen$set_uniform_parameter("q_thresh", lower = 0.90, upper = 0.99, decimals = 2) 

# Growth rate
lhs_gen$set_uniform_parameter("rmax_quantile", lower = 0, upper = 1, decimals = 2)
lhs_gen$set_uniform_parameter("standard_deviation", lower = 0.00, upper = 0.70, decimals = 2)

# Dispersal
lhs_gen$set_uniform_parameter("dispersal_p", lower = 0.05, upper = 0.40, decimals = 2)
## mean dispersal between 5 and 40 km
lhs_gen$set_uniform_parameter("dispersal_b", lower = 5, upper = 40, decimals = 0)
lhs_gen$set_uniform_parameter("dispersal_n_k_threshold", lower = 0.7, upper = 1.0, decimals = 2)

# Density max
## Density: animals/km2 needs to be scaled by grid size (10km x 10km)
## e.g. 1/km2 = (1 animal/km2 * (10*10) ) * frac_cell_used
## 1 km2 = 80 per grid cell = (1*(10*10))*0.8 # assuming 80% grid cell used
## Here I have assumed only 80% of cell is suitable. Upper/lower = 1/km - 6.25/km
lhs_gen$set_uniform_parameter("max_dens", lower = 80, upper = 500, decimals = 0)

# Translocation
lhs_gen$set_uniform_parameter("trans_n", lower = 10, upper = 100, decimals = 0)

sample_data <- lhs_gen$generate_samples(number = 10, random_seed = 42)
head(sample_data)

# Make unique row names for saving files
{set.seed(54612)
  sample_data$UniqueID <- paste0(
    stri_rand_strings(nrow(sample_data), 4, "[A-Z]"),
    stri_rand_strings(nrow(sample_data), 4, "[0-9]"))
}
sample_data <- sample_data[, c(9, 1:8)]
sample_data

We we can run a simulation for each set (or row) of sampled parameters. The SimulationManager class manages the generation of parameters (via the generators), the running the model simulations, and writing simulation results to disk. It also maintains a log of each simulation's success and any errors or warnings encountered.

model <- model_template$clone()
model$set_attributes(params = list("standard_deviation" = NULL,
                                   "dispersal_source_n_k$threshold" = NULL,
                                   "dispersal_source_n_k$cutoff" = 0.00))

# Build the simulation manager
sim_manager <- SimulationManager$new(
  sample_data = sample_data,
  model_template = model,
  # initial_hs = initial_hs,
  generators = list(dispersal_gen, capacity_gen, rmax_gen),
    parallel_cores = 1L,
  results_filename_attributes =
    c(NULL, "UniqueID", "results"),
  results_ext = ".RDS",
  results_dir = OUTPUT_DIR)

# Takes <10 seconds to run 10 example sims on a single core.
system.time({
  run_output <- sim_manager$run()
})
run_output$summary

Step 6: Extract results from simulations

Now that the simulations have run, we can extract the modelled abundances. We now wish to collate summary results for each of our simulations via the ResultsManager class. This manager loads the results from each sample simulation into an intermediate PopulationResults class object, which dynamically generates further results. We need to define functions for calculating summary metrics, as well as any matrices (one row of values per simulation) that we may be interested in examining.

We can see from above that the run_output$summary shows that there were a number of simulations that didn't complete successfully:

run_output$summary

The run_output$full_log shows that some simulations produced NA values when calculating population abundances:

"Warning: Non-finite stage abundances returned by user-defined translocation function"

These simulations had a range of parameter values that caused the simulation to "fall-over". These errors could be fixed by improving the translocation function, or alternatively the user could discard the simulations as being structurally wrong - i.e. that specific combination of parameters is simply unsuitable. We're going to treat them as the latter.

Extract results

# Extract timeseries of abundance from each of the sims
# Load our results (list) into a PopulationResults object
p_results <- PopulationResults$new(results = run_output)
res_manager <- ResultsManager$new(simulation_manager = sim_manager,
                                  simulation_results = p_results,
                                  generators = NULL,
                                  summary_matrices = c("n",
                                                       "distr_pop"),
                                  summary_functions = list(
                                    # total pop abundance
                                    "n" = function(sim_results) {
                                      sim_results$all$abundance
                                    },
                                    # matrix of abundance
                                    ## can be made into raster
                                    "distr_pop" = function(sim_results) {
                                     sim_results$abundance
                                    }
                                  ),
                                  parallel_cores = 1L)
gen_log <- res_manager$generate()
gen_log$summary

# matrix of total population abundances
## each row is a sim, each column a timestep
res_manager$summary_matrix_list$n

# plot
matplot(x = 1:ncol(res_manager$summary_matrix_list$n), 
        y = t(res_manager$summary_matrix_list$n), type = "b",
        lty = 1, xlab = "timestep", ylab = "total abundance")

From the results we can see that our translocation function was working correctly, and that our capacity_gen was defined correctly, as the abundances in the first time step are the same as the trans_n values from our latin-hypercube samples.

identical(unlist(res_manager$summary_matrix_list$n[, 1]),
          unlist(sample_data$trans_n))

Let's assume that samples 2, 3, 4, and 8 are our "best" simulations after some sort of validation. We have already extracted the abundances from these simulations, so now with a few more lines of code, we can generate averages of the simulations as rasters

best_sims <- c(2:4, 8)
dim(res_manager$summary_matrix_list$distr_pop[best_sims, ])

best_abund <- matrix(
  nrow = region$region_cells,
  ncol = 11, # 11 timesteps,
  data = round(colMeans(res_manager$summary_matrix_list$distr_pop[best_sims, ])))

best_abund <- region$raster_from_values(best_abund)
best_abund[[c(1, 6, 11)]]

abd_max <- round_any(max(values(best_abund), na.rm = TRUE),
                     accuracy = 100, ceiling)

# plot of log(x+1) abundances
plot(log1p(best_abund),
     col = hcl.colors(100, "Spectral", rev = TRUE),
     colNA = "grey75",
     addfun = function() {
       plot(ibra_polygons, border = "#000000", col = NA, add = TRUE)
       },
     zlim = c(0, log1p(abd_max)))

poems documentation built on Oct. 7, 2023, 9:06 a.m.