Different types of sensitivity analyses can be conducted using the
nlrx package. To perform a local sensitivity analysis, we recommend
using the simdesign_distinct() to specify local changes of
parameters. Afterwards the proportion of output change can be easily
calculated from the simulation results. The nlrx package also provides
simdesign helper functions to conduct more sophisticated methods such as
Morris Elementary Effects Screening (simdesign_morris()),
Sobol variance decomposition (simdesign_sobol(),
simdesign_sobol2007(),
simdesign_soboljansen()) and Extended Fourier amplitude
sensitivity test (simdesign_eFAST). Additionally, output of
Latin Hypercube Sampling designs (simdesign_lhs()) can be
used to calculate parameter effects based on Partial (rank) correlation
coefficients or Standardised (rank) regression coefficients.
In this vignette, we present an example of the Morris Elementary Effects screening. Other sensitivity analyses simdesigns work in a quite similar way. Details on the specific methods can be found in the corresponding simdesign help pages and the documentation of the sensitivity package. The second example shows how Latin Hypercube Sampling can be used to calculate Partial (rank) correlation coefficients and Standardised (rank) regression coefficients.
Here we present a simple example for running a Morris Sensitivity Analysis with nlrx. We use the Wolf Sheep Predation model from the models library for this example.
library(nlrx)
# Windows default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("C:/Program Files/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("C:/out")
# Unix default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("/home/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("/home/out")
nl <- nl(nlversion = "6.0.3",
nlpath = netlogopath,
modelpath = modelpath,
jvmmem = 1024)In this example, we want to calculate sensitivity for 3 outputs
(number of sheep, number of wolves, number of grass patches). We vary
all numeric model parameters to estimate their sensitivity on the three
defined output metrics. Thus, we define parameter ranges and
distribution functions for all our numeric model parameters. We set the
runtime of the model to 500 ticks and measure our metrics on each tick
(tickmetrics = "true"). However, for calculation of
sensitivity indices, we only want to consider the last 200 ticks. Thus,
we set evalticks to seq(300,500).
nl@experiment <- experiment(expname = "wolf-sheep-morris",
outpath = outpath,
repetition = 1,
tickmetrics = "true",
idsetup = "setup",
idgo = "go",
runtime = 500,
evalticks = seq(300,500),
metrics=c("count sheep", "count wolves", "count patches with [pcolor = green]"),
variables = list("initial-number-sheep" = list(min=50, max=150, step=10, qfun="qunif"),
"initial-number-wolves" = list(min=50, max=150, step=10, qfun="qunif"),
"grass-regrowth-time" = list(min=0, max=100, step=10, qfun="qunif"),
"sheep-gain-from-food" = list(min=0, max=50, step=10, qfun="qunif"),
"wolf-gain-from-food" = list(min=0, max=100, step=10, qfun="qunif"),
"sheep-reproduce" = list(min=0, max=20, step=5, qfun="qunif"),
"wolf-reproduce" = list(min=0, max=20, step=5, qfun="qunif")),
constants = list("model-version" = "\"sheep-wolves-grass\"",
"show-energy?" = "false"))We use the simdesgin_morris() function to attach a
Morris Sensitivity Analysis design. The morrislevels
parameter sets the number of different values for each parameter
(sampling density). The morrisr paramater sets the number
of repeated samplings (sampling size). The morrisgridjump
parameter sets the number of levels that are increased/decreased for
computing the elementary effects. Morris recommendation is to set this
value to levels / 2. We can increase the
nseeds parameter in order to perform multiple runs of the
same parameter matrix with different random seeds. The variation between
those repetitions is an indicator of the stochasticity effects within
the model. More information on the Morris specific parameters can be
found in the description of the morris function in the sensitivity
package (?morris).
nl@simdesign <- simdesign_morris(nl=nl,
morristype="oat",
morrislevels=4,
morrisr=1000,
morrisgridjump=2,
nseeds=5)To execute the simulations, we can use the function
run_nl_all(). Sensitivity analyses typically have many runs
that need to be simulated, thus we recommend to parallelize model runs
by adjusting the future plan (more details on parallelization can be
found in the “Advanced configuration” vignette).
library(future)
plan(multisession)
results <- run_nl_all(nl)First, we need to attach the results to the nl object.
setsim(nl, "simoutput") <- results
saveRDS(nl, file.path(nl@experiment@outpath, "morris.rds"))After results have been attached, we can use the
analyze_nl() function to calculate morris sensetivity
indices.
morris <- analyze_nl(nl)Here we perform a Latin Hypercube Sampling to calculate Partial (rank) correlation coefficients and Standardised (rank) regression coefficients.
library(nlrx)
# Windows default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("C:/Program Files/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("C:/out")
# Unix default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("/home/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("/home/out")
nl <- nl(nlversion = "6.0.3",
nlpath = netlogopath,
modelpath = modelpath,
jvmmem = 1024)In this example, we want to calculate sensitivity for 3 outputs
(number of sheep, number of wolves, number of grass patches). We vary
all numeric model parameters to estimate their sensitivity on the three
defined output metrics. Thus, we define parameter ranges and
distribution functions for all our numeric model parameters. We set the
runtime of the model to 500 ticks and measure our metrics on each tick
(evalticks = "true").
nl@experiment <- experiment(expname = "wolf-sheep-morris",
outpath = outpath,
repetition = 1,
tickmetrics = "true",
idsetup = "setup",
idgo = "go",
runtime = 500,
metrics=c("count sheep", "count wolves", "count patches with [pcolor = green]"),
variables = list("initial-number-sheep" = list(min=50, max=150, step=10, qfun="qunif"),
"initial-number-wolves" = list(min=50, max=150, step=10, qfun="qunif"),
"grass-regrowth-time" = list(min=0, max=100, step=10, qfun="qunif"),
"sheep-gain-from-food" = list(min=0, max=50, step=10, qfun="qunif"),
"wolf-gain-from-food" = list(min=0, max=100, step=10, qfun="qunif"),
"sheep-reproduce" = list(min=0, max=20, step=5, qfun="qunif"),
"wolf-reproduce" = list(min=0, max=20, step=5, qfun="qunif")),
constants = list("model-version" = "\"sheep-wolves-grass\"",
"show-energy?" = "false"))Here we want to run a Latin Hypercube Sampling, thus we use the
simdesign_lhs() function.
nl@simdesign <- simdesign_lhs(nl, samples=500, nseeds=1, precision=3)To execute the simulations, we can use the function
run_nl_all(). Sensitivity analyses typically have many runs
that need to be simulated, thus we recommend to parallelize model runs
by adjusting the future plan (more details on parallelization can be
found in the “Advanced configuration” vignette).
library(future)
plan(multisession)
results <- run_nl_all(nl, split=10)First, we need to attach the results to the nl object.
setsim(nl, "simoutput") <- results
saveRDS(nl, file.path(nl@experiment@outpath, "lhs.rds"))After results have been attached, we need to post-process our data to
run the pcc and src function of the
sensitivity package. We first take our parameter matrix
(siminput) and select only columns with variable parameters
and drop all other columns. We also need to rename the columns because
pcc and src do not support special characters
(-) in column names.
Our simulation results are measured for each tick, thus we first need to aggregate our output. Here we just calculate the mean and standard deviation of outputs for each random-seed and siminputrow combination. Afterwards, we drop the random seed and siminputrow columns and rename the columns to remove special characters.
Finally, we use both datasets to run the pcc and
src functions. These functions can only compute
coefficients for one output at a time. Thus, we nested the function call
inside a purrr::map() function that iterates over the
column names of our output tibble.
library(tidyverse)
input <- getsim(nl, "siminput") %>% # Take input parameter matrix
dplyr::select(names(getexp(nl, "variables"))) %>% # Select variable parameters only
dplyr::rename_all(~str_replace_all(., c("-" = "_", "\\s+" = "_"))) # Remove - and space characters.
output <- getsim(nl, "simoutput") %>% # Take simulation output
dplyr::group_by(`random-seed`, siminputrow) %>% # Group by random seed and siminputrow
dplyr::summarise_at(getexp(nl, "metrics"), list(mean=mean, sd=sd)) %>% # Aggregate output
dplyr::ungroup() %>% # Ungroup
dplyr::select(-`random-seed`, -siminputrow) %>% # Only select metrics
dplyr::rename_all(~str_replace_all(., c("-" = "_", "\\s+" = "_", "\\[" = "_", "\\]" = "_", "=" = ""))) # Remove - and space characters.
# Perform pcc and src for each output separately (map)
pcc.result <- purrr::map(names(output), function(x) sensitivity::pcc(X=input, y=output[,x], nboot = 100, rank = FALSE))
src.result <- purrr::map(names(output), function(x) sensitivity::src(X=input, y=output[,x], nboot = 100, rank = FALSE)) The results are reported as a nested list, where each outer element
represents one of the calculated model outputs. The inner list items
represent the different outputs from the pcc and
src functions.
We can for example look at the pcc results of one
specific output by using the basic plot function:
plot(pcc.result[[1]])We can also extract all the data to a tidy data format and create nice plots with the ggplot package:
pcc.result.tidy <- purrr::map_dfr(seq_along(pcc.result), function(x) {
pcc.result[[x]]$PCC %>%
tibble::rownames_to_column(var="parameter") %>%
dplyr::mutate(metric = names(output)[x])
})
ggplot(pcc.result.tidy) +
coord_flip() +
facet_wrap(~metric) +
geom_point(aes(x=parameter, y=original, color=metric)) +
geom_errorbar(aes(x=parameter, ymin=`min. c.i.`, ymax=`max. c.i.`, color=metric), width=0.1)
src.result.tidy <- purrr::map_dfr(seq_along(src.result), function(x) {
src.result[[x]]$SRC %>%
tibble::rownames_to_column(var="parameter") %>%
dplyr::mutate(metric = names(output)[x])
})
ggplot(src.result.tidy) +
coord_flip() +
facet_wrap(~metric) +
geom_point(aes(x=parameter, y=original, color=metric)) +
geom_errorbar(aes(x=parameter, ymin=`min. c.i.`, ymax=`max. c.i.`, color=metric), width=0.1)