scenarios: Compare epidemic scenarios

scenarios is intended to provide functions to compare the outcomes of epidemic modelling simulations. This package is still a work in progress.

Installation

You can install the development version of scenarios from GitHub with:

if (!require("remotes")) install.packages("remotes")
remotes::install_github("epiverse-trace/scenarios")

Quick start

An example with finalsize

Define an epidemic model scenario by creating a new scenario object. The standard workflow needs a model function, such as finalsize::final_size(), and appropriate arguments to the function.

# load scenarios
library(scenarios)

# create a scenario for pandemic-potential influenza
# with finalsize::final_size() as the model function
scenario_pandemic_flu <- scenario(
  model_function = "finalsize::final_size",
  parameters = make_parameters_finalsize_UK(), # a helper function
  replicates = 3
)

View a summary of the scenario.

scenario_pandemic_flu
#> Epidemic scenario object
#>  Model function: finalsize::final_size
#>  Scenario replicates: 3
#>  Scenario outcomes are not prepared

The scenario object is created but the model function is not initially run. This can be checked using a helper function. Tip: Many helper functions have the prefix sce_ to help them be found quickly when using autocomplete in various text editors and IDEs.

# check whether the scenario data have been generated
sce_has_data(scenario_pandemic_flu)
#> [1] FALSE

Model outcome data can be generated by running the model function with the parameters specified in the scenario.

scenario_pandemic_flu <- run_scenario(scenario_pandemic_flu)

Take a peek at the column names in the model outcome replicates.

sce_peek_outcomes(scenario_pandemic_flu)
#>       demo_grp       susc_grp susceptibility     p_infected      replicate 
#>    "character"    "character"      "numeric"      "numeric"      "integer"

Get the outcomes from all replicates as a single dataset, or aggregate an outcome variable of interest across replicates by some grouping variable.

# get all output
head(sce_get_outcomes(scenario_pandemic_flu))
#>    demo_grp   susc_grp susceptibility p_infected replicate
#> 1:   [0,20) susc_grp_1              1  0.6544866         1
#> 2:  [20,40) susc_grp_1              1  0.5750030         1
#> 3:      40+ susc_grp_1              1  0.4588871         1
#> 4:   [0,20) susc_grp_1              1  0.6544866         2
#> 5:  [20,40) susc_grp_1              1  0.5750030         2
#> 6:      40+ susc_grp_1              1  0.4588871         2

# aggregate proportion infected by demographic group
# NOTE that all replicates have the same outcome in this deterministic model
sce_aggregate_outcomes(
  x = scenario_pandemic_flu,
  grouping_variables = c("demo_grp"),
  measure_variables = c("p_infected"),
  summary_functions = c("mean", "min", "max")
)
#>    demo_grp p_infected_mean p_infected_min p_infected_max
#> 1:      40+       0.4588871      0.4588871      0.4588871
#> 2:   [0,20)       0.6544866      0.6544866      0.6544866
#> 3:  [20,40)       0.5750030      0.5750030      0.5750030

An example with epidemics

This example shows the same workflow applied to a simple, deterministic epidemic model from the epidemics package.

# create a new scenario
scenario_sir <- scenario(
  model_function = "epidemics::sir_desolve",
  parameters = make_parameters_SIR_epidemic(), # a helper function
  replicates = 5L
)

# view the initial conditions and infection parameters
sce_get_information(scenario_sir, which = c("init", "parms"))
#> $init
#>    S    I    R 
#> 0.99 0.01 0.00 
#> 
#> $parms
#>  beta gamma 
#>   1.0   0.1

# generate scenario outcomes by running the model
scenario_sir <- run_scenario(scenario_sir)

# peek at the outcomes
sce_peek_outcomes(scenario_sir)
#>       time      state proportion  replicate 
#>  "numeric"   "factor"  "numeric"  "integer"

# view the aggregated outcomes
# this is the per-timestep, per-class (S, I, R) mean proportion
# and is the same across replicates in this deterministic model
tail(
  sce_aggregate_outcomes(
    scenario_sir,
    grouping_variables = c("time", "state"),
    measure_variables = "proportion",
    summary_functions = "mean"
  )
)
#>    time state proportion_mean
#> 1:   99     S    4.501519e-05
#> 2:   99     I    8.643038e-05
#> 3:   99     R    9.998686e-01
#> 4:  100     S    4.501149e-05
#> 5:  100     I    7.820928e-05
#> 6:  100     R    9.998768e-01

Help

To report a bug please open an issue.

Contribute

Contributions to scenarios are welcomed. Please follow the package contributing guide.

Code of conduct

Please note that the scenarios project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.