inst/doc/sensitivity.R

In nlrx: Setup, Run and Analyze 'NetLogo' Model Simulations from 'R' via 'XML'

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
#  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)

## ----eval=FALSE---------------------------------------------------------------
#  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"))

## ----eval=FALSE---------------------------------------------------------------
#  nl@simdesign <- simdesign_morris(nl=nl,
#                                   morristype="oat",
#                                   morrislevels=4,
#                                   morrisr=1000,
#                                   morrisgridjump=2,
#                                   nseeds=5)

## ----eval=FALSE---------------------------------------------------------------
#  library(future)
#  plan(multisession)
#  results <- run_nl_all(nl)

## ----eval=FALSE---------------------------------------------------------------
#  setsim(nl, "simoutput") <- results
#  saveRDS(nl, file.path(nl@experiment@outpath, "morris.rds"))

## ----eval=FALSE---------------------------------------------------------------
#  morris <- analyze_nl(nl)

## ----eval=FALSE---------------------------------------------------------------
#  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)

## ----eval=FALSE---------------------------------------------------------------
#  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"))

## ----eval=FALSE---------------------------------------------------------------
#  nl@simdesign <- simdesign_lhs(nl, samples=500, nseeds=1, precision=3)

## ----eval=FALSE---------------------------------------------------------------
#  library(future)
#  plan(multisession)
#  results <- run_nl_all(nl, split=10)

## ----eval=FALSE---------------------------------------------------------------
#  setsim(nl, "simoutput") <- results
#  saveRDS(nl, file.path(nl@experiment@outpath, "lhs.rds"))

## ----eval=FALSE---------------------------------------------------------------
#  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))

## ----eval=FALSE---------------------------------------------------------------
#  plot(pcc.result[[1]])

## ----eval=FALSE---------------------------------------------------------------
#  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)