inst/doc/abc.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",
#                              outpath=outpath,
#                              repetition=1,
#                              tickmetrics="false",
#                              idsetup="setup",
#                              idgo="go",
#                              runtime=100,
#                              metrics=c("count sheep", "count wolves"),
#                              variables = list("sheep-gain-from-food" = list(min=2, max=6, qfun="qunif"),
#                                               "wolf-gain-from-food" = list(min=10, max=30, qfun="qunif")),
#                              constants = list('initial-number-sheep' = 100,
#                                               'initial-number-wolves' = 50,
#                                               "grass-regrowth-time" = 30,
#                                               "sheep-reproduce" = 4,
#                                               "wolf-reproduce" = 5,
#                                               "model-version" = "\"sheep-wolves-grass\"",
#                                               "show-energy?" = "false"))
#  

## ----eval=FALSE---------------------------------------------------------------
#  nl@simdesign <- simdesign_ABCmcmc_Marjoram(nl=nl,
#                                             summary_stat_target = c(100, 80),
#                                             n_rec = 100,
#                                             n_calibration=200,
#                                             use_seed = TRUE,
#                                             progress_bar = TRUE,
#                                             nseeds = 1)

## ----eval=FALSE---------------------------------------------------------------
#  post <- function(nl){
#    res <- getsim(nl, "simoutput") %>%
#      dplyr::select(getexp(nl, "metrics")) %>%
#      dplyr::summarise_each(list(max=max))
#    return(as.numeric(res))
#  }
#  
#  nl@simdesign <- simdesign_ABCmcmc_Marjoram(nl=nl,
#                                             postpro_function = post,
#                                             summary_stat_target = c(100, 80),
#                                             n_rec = 100,
#                                             n_calibration=200,
#                                             use_seed = TRUE,
#                                             progress_bar = TRUE,
#                                             nseeds = 1)

## ----eval=FALSE---------------------------------------------------------------
#  results <- run_nl_dyn(nl, seed = nl@simdesign@simseeds[1])

## ----eval=FALSE---------------------------------------------------------------
#  setsim(nl, "simoutput") <- results
#  saveRDS(nl, file.path(nl@experiment@outpath, "ABCmcmc.rds"))
#  
#  ## Calculate descriptive statistics
#  getsim(nl, "simoutput") %>% # get simulation results from nl object
#    dplyr::select(param) %>% # select param column
#    tidyr::unnest(cols=param) %>%  # unnest param column
#    dplyr::summarise_each(list(min=min, max=max, mean=mean, median=median)) %>% # calculate statistics
#    tidyr::gather(parameter, value) %>% # convert to long format
#    tidyr::separate(parameter, into = c("parameter", "stat"), sep = "_") %>% # seperate parameter name and statistic
#    tidyr::spread(stat, value) # convert back to wide format
#  
#  ## Plot histogram of parameter sampling distribution:
#  getsim(nl, "simoutput") %>% # get simulation results from nl object
#    dplyr::select(param) %>% # select param column
#    tidyr::unnest(cols=param) %>%  # unnest param column
#    tidyr::gather(parameter, value) %>% # convert to long format
#    ggplot2::ggplot() + # plot histogram with a facet for each parameter
#    ggplot2::facet_wrap(~parameter, scales="free") +
#    ggplot2::geom_histogram(ggplot2::aes(x=value), bins = 40)
#  
#  ## Plot density of parameter sampling distribution:
#  getsim(nl, "simoutput") %>% # get simulation results from nl object
#    dplyr::select(param) %>% # select param column
#    tidyr::unnest(cols=param) %>% # unnest param column
#    tidyr::gather(parameter, value) %>% # convert to long format
#    ggplot2::ggplot() + # plot density with a facet for each parameter
#    ggplot2::facet_wrap(~parameter, scales="free") +
#    ggplot2::geom_density(ggplot2::aes(x=value, fill=parameter))
#  
#  ## Get best parameter combinations and corresponding function values
#  getsim(nl, "simoutput") %>%  # get simulation results from nl object
#    dplyr::select(dist,epsilon) %>%  # select dist and epsilon columns
#    tidyr::unnest(cols=c(dist,epsilon)) %>%  # unnest dist and epsilon columns
#    dplyr::mutate(runID=dplyr::row_number()) %>% # add row ID column
#    dplyr::filter(dist == epsilon) %>% # only keep runs with dist=epsilon
#    dplyr::left_join(getsim(nl, "simoutput") %>% # join parameter values of best runs
#                       dplyr::select(param) %>%
#                       tidyr::unnest(cols=param) %>%
#                       dplyr::mutate(runID=dplyr::row_number())) %>%
#    dplyr::left_join(getsim(nl, "simoutput") %>% # join output values best runs
#                       dplyr::select(stats) %>%
#                       tidyr::unnest(cols=stats) %>%
#                       dplyr::mutate(runID=dplyr::row_number())) %>%
#    dplyr::select(runID, dist, epsilon, dplyr::everything()) # update order of columns
#  
#  ## Analyse mcmc using coda summary and plot functions:
#  summary(coda::as.mcmc(getsim(nl, "simoutput") %>%
#                          dplyr::select(param) %>%
#                          tidyr::unnest(cols=param)), quantiles =c(0.05,0.95,0.5))
#  
#  plot(coda::as.mcmc(getsim(nl, "simoutput") %>%
#                          dplyr::select(param) %>%
#                          tidyr::unnest(cols=param)))
#  

## ----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",
#                              outpath=outpath,
#                              repetition=1,
#                              tickmetrics="false",
#                              idsetup="setup",
#                              idgo="go",
#                              runtime=100,
#                              metrics=c("count sheep", "count wolves"),
#                              variables = list("sheep-gain-from-food" = list(min=2, max=6, qfun="qunif"),
#                                               "wolf-gain-from-food" = list(min=10, max=30, qfun="qunif")),
#                              constants = list('initial-number-sheep' = 100,
#                                               'initial-number-wolves' = 50,
#                                               "grass-regrowth-time" = 30,
#                                               "sheep-reproduce" = 4,
#                                               "wolf-reproduce" = 5,
#                                               "model-version" = "\"sheep-wolves-grass\"",
#                                               "show-energy?" = "false"))
#  

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

## ----eval=FALSE---------------------------------------------------------------
#  results <- run_nl_all(nl)

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

## ----eval=FALSE---------------------------------------------------------------
#  input <- getsim(nl, "siminput") %>%
#    dplyr::select(names(getexp(nl, "variables")))
#  output <- getsim(nl, "simoutput") %>%
#    dplyr::select(getexp(nl, "metrics"))
#  target <- c("count sheep"=100, "count wolves"=80)

## ----eval=FALSE---------------------------------------------------------------
#  results.abc.reject <- abc::abc(target=target,
#                          param=input,
#                          sumstat=output,
#                          tol=0.3,
#                          method="rejection")
#  
#  results.abc.loclin <- abc::abc(target=target,
#                                 param=input,
#                                 sumstat=output,
#                                 tol=0.3,
#                                 method="loclinear")

## ----eval=FALSE---------------------------------------------------------------
#  results.abc.all <- tibble::as_tibble(results.abc.reject$unadj.values) %>% # results from rejection method
#    tidyr::gather(parameter, value) %>%
#    dplyr::mutate(method="rejection") %>%
#    dplyr::bind_rows(tibble::as_tibble(results.abc.loclin$adj.values) %>% # results from local linear regression method
#                       tidyr::gather(parameter, value) %>%
#                       dplyr::mutate(method="loclinear")) %>%
#    dplyr::bind_rows(input %>%                # initial parameter distribution (lhs)
#                       tidyr::gather(parameter, value) %>%
#                       dplyr::mutate(method="lhs"))
#  
#  ggplot2::ggplot(results.abc.all) +
#    ggplot2::facet_grid(method~parameter, scales="free") +
#    ggplot2::geom_histogram(ggplot2::aes(x=value))
#  
#  ggplot2::ggplot(results.abc.all) +
#    ggplot2::facet_wrap(~parameter, scales="free") +
#    ggplot2::geom_density(ggplot2::aes(x=value, fill=method), alpha=0.1)
#