Models are nowadays so important and widespread that we often don’t realize that most of what we use in our daily lives is the result of someone’s (or some team’s) model. From weather forecasting, to predicting stock market fluctuations. From the design of your bike, to the estimation of blood sugar and insulin levels and how they vary hourly and according to different diets (in case you don’t know, this can be important for diabetes treatments, and drug dosage and scheduling)… Models are indeed ‘everywhere’.
In a world where disciplines often interact to study complex questions, like climate change, natural resource management, or even where to build an offshore wind farm, modellers (and model users!) face challenges that can impair scientific and social progress:
Continuous adaptive management is an existing framework that assumes a periodic re-assessment of the status of a situation or issue, e.g. the amount of waste a company produces. This leads to the design and implementation of new solutions - e.g. recycling paper or using rain water - their monitoring and evaluation, and their adjustment if need be, before the next assessment.
In modelling terms, this framework implies that for a given problem/question analyses and forecasts are repeated as new data, new methodologies and new or improved models arise. This can lead to a re-assessment of risk, inform decisions and influence policy making. These can in turn feedback into the process of re-framing the initial problem/question.
Although this framework is ideal, it is only seldomly applied across the scientific community.
Why? Because we are yet far from developing our research in a way that is easily repeatable, “upgradeable” and shareable. This renders the application of the continuous adaptive management framework costly, both temporally and financially, especially in interdisciplinary projects that involve multiple teams and stakeholders.
SpaDES - Spatial Discrete Event Simulator - is a modelling framework in the form of an ensemble of
R packages. It aims to provide modellers and non-modellers across different domains of life and earth sciences (and beyond!) with a modelling tool that they can control, apply and develop. In its whole,
SpaDES bridges the gaps between modellers and model-users of different disciplines, between scientists, practitioners, managers, stakeholders and policy makers, and between scientists and the general public.
SpaDES is written in
R: a widely understood programming language, highly suited to statistical analyses, but flexible and fast enough to accommodate simulation work;
SpaDES can be used for a wide-range of modelling approaches: + process-based models; + phenomenological models; + spatially-explicit models; + … 3.
SpaDES models are built of modules, i.e. self-contained, interacting pieces that execute a particular event, or process in the model; 4.
SpaDES modules can be made available across the entire
R community if module-developers wish it so; 5.
SpaDES will soon include a web-app builder to make interactive simulation possible.
Like a smart-phone where users install and use the apps that they want,
SpaDES can be seen as a platform where a model can be “assembled” from various pre-made modules. This allows non-modellers to run someone else’s model on their data and their parameters, and, on the other hand, it allows modellers write their own modules from scratch, run them and share them with the
R user community.
Our question(s): how will forest fires affect forest age? And how will climate change affect fire regimes and in turn the vegetation? How will timber supply be affected?
To start we’ll need:
* a vegetation dynamics module - e.g. forestAge; * a fire simulator - e.g. forestSuccessionBeacons modules; * optional: + the interaction between fire and vegetation + statistical analyses/presentation of results
Depending on how modules are built, the interaction between modules may or may not require additional modules. Similarly, the presentation of results and their analyses can be directly made within the module producing them, or as posteriori using another module. Your call.
SpaDES is all about flexibility!
Using the above mentioned modules would already allow an analysis of how fires affect forest age, and vice-versa.
Then add climate-change, timber supply models…
Analysing the effects of climate change could then be done by varying fire-related parameters (e.g. frequency and intensity of fires) in a simulation experiment. Alternatively, an additional module could simulate climate change effects, which would feedback unto both fire parameters and vegetation.
Finally, a timber supply module could “read” the vegetation maps produced after fire simulations and calculate timber quantity and quality from forest conditions (e.g. surviving tree species and their age).
Reserve, add new data, and re-run
Provided that modules are able to download, read and treat data adequately, updating the simulations for new data is easy.
SpaDES comes with a set of functions (via the
reproducible package) that enable caching model structure, parameters and simulation results. This way, if new data only affects a section of the model, the user won’t have to repeat every step to update the simulations.
* Learn more about: sharing modules and caching (here and here)
Decorate and serve
Because models are not just for modellers, a shiny app building package is being developed to transform any
SpaDES simulation into a web-app. This provides an interactive way of visualising results, as users can change simulation parameters and turn modules on and off, whilst observing how this impacts simulation results. Importantly, it can be particularly useful when conveying results to a non-expert or even the public.
SpaDES groups a series of
R packages - metapackage - containing both core and additional methods to design a discrete event simulation model:
R- having R Studio, or another
Reditor might also be useful.
SpaDESis still on a relatively steep development curve. As such, it might be sometimes necessary to update a particular package from Git Hub (where
SpaDES is being developed). It might be important to keep an eye on new improvements present in Git Hub that might not be in CRAN yet. You’ll find the latest releases for each package in their respective Git Hub page (e.g. SpaDES.core releases). If you go to the commits page you’ll notice that that there might be new features relatively to the latest commit.
## Install all packages and their dependencies first install.packages('SpaDES', dependencies = TRUE) ## now, update some of the packages from the development branch. devtools::install_github("PredictiveEcology/quickPlot", ref = "CEFWorkshop" , dependencies = FALSE) devtools::install_github("PredictiveEcology/reproducible", ref = "CEFWorkshop", dependencies = FALSE) devtools::install_github("PredictiveEcology/SpaDES.tools", ref = "CEFWorkshop", dependencies = FALSE) devtools::install_github("PredictiveEcology/SpaDES.core", ref = "CEFWorkshop", dependencies = FALSE)
## load SpaDES library(SpaDES) ## Locate the sample modules in SpaDES.core sampleModuleDir <- system.file('sampleModules', package = 'SpaDES.core') ## list the modules in this directory dir(sampleModuleDir)
What do these modules do?
## open the sample module Rmd file to learn more about the modules and access a simple example sampleModuleFiles <- list.files(sampleModuleDir, pattern = 'Rmd$', recursive = TRUE, full.names = TRUE) ## open the first file file.edit(sampleModuleFiles) ## or all: lapply(sampleModuleFiles, file.edit)
All modules should come with an .Rmd file that provides a summary of how the module works and reproducible example. If you’ve ever searched for R-help on stackoverflow you’ll have noticed how people try to provide
R code that can be run in any machine to reproduce the problem/solution at hand. Module .Rmd files should also run in any machine.
Let’s go back to the modelling question discussed briefly above.
Our first objective is to simulate the effect of fire on forest age. There are two main components to such a model, which need to interact:
LCC2005 module group provides a set of modules that can reproduce these dynamics (in a simple way).
Let’s define the working directories of our simulation exercise (SpaDESintro) first.
main.dir <- "~/temp_SpaDESws/" ## We'll create the necessary directories in a temporary location - feel free to change this if you wish to. setPaths(cachePath = file.path(main.dir,'SpaDESintro/cache'), inputPath = file.path(main.dir, 'SpaDESintro/inputs'), modulePath = file.path(main.dir, 'SpaDESintro/modules'), outputPath = file.path(main.dir, 'SpaDESintro/outputs')) ## verify that this is want you wanted: getPaths()
Now we’ll get
LCC2005 modules needed and save them in the
modules folder created before. If you already have done this before and want to do it again, set the
overwrite argument to
data = TRUE will make sure the necessary data for the module is downloaded. Unfortunately, the forest age raster (“can_age04_1km.tif”) is no longer downloadable without creating an account on the NASA EarthData website. Download it from here and put in inside modules/forestAges/data (this will only be possible during the workshop).
For future reference, it is possible to create an account on the NASA EarthData website and, after logging in, you’ll find the file here.
## note how we use the modules folder defined above and chose to download the necessary data downloadModule('LCC2005', path = getPaths()$modulePath, repo = 'CeresBarros/SpaDES-modules', data = TRUE, overwrite = TRUE) ## open the description of the LCC2005 module group - on Windows this provides files directories openModules('LCC2005', path = getPaths()$modulePath) ## list all module files within LCC2005 and open them LCC2005ModuleFiles <- list.files(getPaths()$modulePath, pattern = 'Rmd$', recursive = TRUE, full.names = TRUE) ## open module Rmd files file.edit(LCC2005ModuleFiles) lapply(LCC2005ModuleFiles, file.edit)
As you can see,
LCC2005 is a module group, meaning that it contains several modules that work together. Find more about module groups here.
In this exercise, we will follow a simplified version of the example contained in the
LCC2005.Rmd file. We will only use the modules necessary to simulate vegetation dynamics (
LccToBeaconsReclassify) and fire spread (
## We've defined the paths earlier paths <- getPaths() ## we'll use several, but not all the modules in LCC2005 modules <- list("cropReprojectLccAge", "fireSpreadLcc", "forestAge", "forestSuccessionBeacons", "LccToBeaconsReclassify") ## we don't need to supply external objects objects <- list() ## define the starting and end times of the simulation times <- list(start = 1.0, end = 20.0) ## define parameters. This will override any default parameters during this simulation. parameters <- list( .globals = list(stackName = "landscape", burnStats = "nPixelsBurned"), .progress = list(NA), fireSpreadLcc = list(nFires = 10L, drought = 1, its = 1e6, persistprob = 0, returnInterval = 1, startTime = 0, .plotInitialTime = NA, .plotInterval = 10) ) ## list parameters ls.str(parameters) ## make simulation object mySim <- simInit(times = times, params = parameters, modules = modules, objects = objects, paths = paths)
mySim contains all the necessary information to run the simulation, such as the objects that are shared between modules, and how these are linked between them. It’s called a
simList. You can find out more about the
simList object here.
Let’s have a look at what mySim contains and how modules are connected.
## list the contents of mySim ls.str(mySim) dev() plot(mySim$age) ## plot module connections dev() moduleDiagram(mySim) objectDiagram(mySim) ## see which events are scheduled events(mySim)
You’ll note that
SpaDES has scheduled the initialisation events for each module.
Now that we’ve checked the structure of our simulation model, we can go ahead and run it.
dev() # to open a new plot device for faster plotting clearPlot() mySim2 <- spades(mySim, debug = TRUE) # set debug = TRUE to print simulation steps. ## After running the simulation, we can print the event diagram, the completed events and the scheduled events eventDiagram(mySim2) completed(mySim2) events(mySim2)
SpaDES it is possible to change a simulations’ parameters and then continue running it. This is possible because, even though the simulation stopped at year 20,
SpaDES has scheduled the “next” events (see the result of the
events call above). Find out more about events and event scheduling here.
What will happen if drought intensity increases and fires become more frequent?
end(mySim2) ## accesses the last year end(mySim2) <- 50 ## add another 30 years to the end of the simulation. ## Alternatively use end(mySim2) <- end(mySim2) + 30 params(mySim2)$fireSpreadLcc$drought ## accesses the current drought parameter value params(mySim2)$fireSpreadLcc$drought <- 3 ## changes it params(mySim2)$fireSpreadLcc$returnInterval ## accesses the current fire return interval value params(mySim2)$fireSpreadLcc$returnInterval <- 1 ## changes it clearPlot() mySim3 <- spades(mySim2, debug = TRUE) time(mySim2) ## "current" time (i.e. the time at which the simulation finished) time(mySim3)
The second part of our question was to investigate how forest fires may affect harvesting. The
LCC2005 module group does not contain any module that we can use for this, but that’s OK, because we can get modules from other Git repositories.
We will use an extremely simple module that maps areas that are adequate for forest harvesting in function of forest age (a parameter that the user inputs) and the type of forest that can be harvested (here, only coniferous forests). This module, called
forestHarvest, will be downloaded from a different repository “CeresBarros/forestHarvest-SpaDESmodule”. Learn more about downloading and sharing modules here.
downloadModule(name = "forestHarvest", repo = "CeresBarros/forestHarvest-SpaDESmodule", path = getPaths()$modulePath, overwrite = TRUE) ## We've defined the paths earlier paths <- getPaths() modules <- list("cropReprojectLccAge", "fireSpreadLcc", "forestAge", "forestSuccessionBeacons", "LccToBeaconsReclassify", "forestHarvest") ## we don't need to suply external objects objects <- list() ## define the starting and end times of the simulation times <- list(start = 1.0, end = 50.0) ## notice the new parameter list parameters <- list( .globals = list(stackName = "landscape", burnStats = "nPixelsBurned"), .progress = list(NA), fireSpreadLcc = list(nFires = 10L, drought = 1, its = 1e6, persistprob = 0, returnInterval = 10, startTime = 0, .plotInitialTime = 0.1, .plotInterval = 10), forestHarvest = list(returnInterval = 5, harvAge = 50, startTime = 11, .plotInterval = 5) ) ## make simulation object mySim_manag <- simInit(times = times, params = parameters, modules = modules, objects = objects, paths = paths) dev() moduleDiagram(mySim_manag) objectDiagram(mySim_manag) dev() clearPlot() mySim_manag2 <- spades(mySim_manag, debug = TRUE)
It is possible to set
SpaDES to save particular objects at a user-defined interval. For this, we need to change the
simInit call to include the
outputs parameter, where we define the objects to be saved, the saving interval, etc. Check the
R help for
simInit for more details.
## Change simInit to include the outputs options mySim_manag <- simInit(times = times, params = parameters, modules = modules, objects = objects, paths = paths, outputs = data.frame(objectName = "ageMap", saveTime = seq(1, 50, by = 10))) ## Re-run the model mySim_manag2 <- spades(mySim_manag, debug = TRUE) ## Print the list of saved outputs outputs(mySim_manag2) ## Read the first object (two) ageMap01 <- readRDS(outputs(mySim_manag2)$file) ## Now you can visualise the saved map in any way you want clearPlot() Plot(ageMap01, title = "Age map") table(ageMap01) Plot(hist(ageMap01), title = "Age distribution")
Are there fewer areas that are appropriate for harvesting if fires become larger (i.e. drought is more severe)?
end(mySim_manag2) <- 100 params(mySim_manag2)$fireSpreadLcc$drought <- 3 dev() clearPlot() mySim_manag3 <- spades(mySim_manag2, debug = TRUE)
SpaDES is still being developed it has already been used to “translate” other important models that have been vastly used in ecology, both for purely academic purposes and in an applied context. One of them is LANDIS II. LANDIS-II is a fairly complex landscape model that simulates forest succession in terms of biomass per species age cohorts (Scheller et al. 2007) and can be coupled to other models to simulate vegetation responses to various drivers, like climate change and fire. What is the interest of translating LANDIS-II into
R? It allows a better control and transparency of all pre-simulation work.
R brings the stats to the model and the work becomes self-contained, transparent and reproducible.
McIntire’s team has already done the hard part - bringing LANDIS into
R - as new users find new applications to the models, they might want to develop other modules that will interact with LBM-R.
## get LBMR module - note: this is temporary as this module is yet unpublished downloadModule("LBMR", repo = "CeresBarros/SpaDES-modules", path = getPaths()$modulePath, data = TRUE, overwrite = TRUE) ## setup the simulation paths <- getPaths() times <- list(start = 0, end = 10) parameters <- list( .globals = list(verbose = FALSE), LBMR = list(.plotInitialTime = 1, successionTimestep = 2) ) modules <- list("LBMR") mySim <- simInit(times = times, params = parameters, modules = modules, paths = paths) dev() objectDiagram(mySim) ## run the simulation clearPlot() mySimOut <- spades(mySim, debug = TRUE)
Like before, parameters can be changed, notably the seed dispersal mechanism. The default is the
wardDispersal, which … But
universalDispersal are also available.
## Check the available parameters parameters(mySim) ## Note how this differs from ls.str(parameters) parameters <- list( .globals = list(verbose = FALSE), LBMR = list(.plotInitialTime = 1, successionTimestep = 2, seedingAlgorithm = "noDispersal") ) modules <- list("LBMR") dev() clearPlot() mySim <- simInit(times = times, params = parameters, modules = modules, paths = paths) mySimOut <- spades(mySim, debug = TRUE)
SpaDES will have another package that will allow building web-apps from
simList objects, so that simulation results can be visualised and changed interactively. For now, you can play around with the proof-of-concept Shiny apps: