Guide to developing epidemics features

Scope

This vignette is structured around use cases of making package modifications given the existing package architecture.

Consequently this vignette does not cover major changes that would require a substantial rewrite of the package, such as switching the C++ solver used for ODE models, or switching to using cpp11 rather than Rcpp.

This vignette also does not tackle questions around epidemic modelling choices - such as which compartmental transitions should be allowed; we strongly recommend the involvement of a specialist in such cases.

Finally, this vignette relates only to the deterministic ODE models in epidemics, as these are expected to be the main starting points for most users seeking to create their own models.

A guide to package structure

A simplified view of the package structure is shown below. Note the files indicated by comments in capitals and the associated ‘level’ of the code; this roughly corresponds to the level in the call stack.

.
├── DESCRIPTION
├── NAMESPACE
├── R
│   ├── RcppExports.R           # auto-generated from `src/`
│   ├── check_args_*.R          # MODEL-SPECIFIC INPUT CROSS-CHECKING FUNC
│   ├── dummy_elements.R        # functions for empty composable elements
│   ├── helpers.R               # functions to handle model output
│   ├── input_check_helpers.R   # functions to construct cross-checking fns
│   ├── intervention.R          # intervention classes
│   ├── model_*.R               # USER-FACING MODEL FUNCTION; LEVEL 1
│   ├── population.R            # population class
│   ├── tools.R                 # internal utility functions
│   └── vaccination.R           # vaccination class
├── epidemics.Rproj
├── inst
│   └── include
│       ├── epidemics.h         # package header, must include all `model_*.h`
│       ├── intervention.h      # functions to handle interventions in C++
│       ├── model_*.h           # MODEL ODE SYSTEM; LEVEL 3
│       ├── ode_tools.h         # definition of the initial state struct
│       ├── population.h        # functions to handle populations in C++
│       ├── time_dependence.h   # functions to handle time-dependence in C++
│       └── vaccination.h       # functions to handle vaccinations in C++
├── man
│   └── *.Rd                    # R function documentation
├── src
│   ├── Makevars                # UNIX make options for Rcpp
│   ├── Makevars.win            # Windows make options for Rcpp
│   ├── RcppExports.cpp         # auto-generated from `src/`
│   └── model_*.cpp             # MODEL FUNC WRAPPING ODE SYSTEM; LEVEL 2
│
├── tests                       # test files
└── vignettes                   # vignettes

Adding or removing model parameters

Adding or removing parameters to the ODE models involves working with all levels of the model code.

Consider an example of adding a uniform, age-independent, background mortality rate ω that is unrelated to the epidemic itself. Removing a parameter involves reversing these steps.

1. Modify the user-facing model function in R/model_*.R to accept the mortality rate as an argument; add documentation and input checks, and ensure that it is included in parameter set creation; this is level 1, the top level, of the code.

   #' @export
   model_new = function(
     population,
     <EXISTING PARAMETERS>,
     mortality_rate,           # <== New parameter added.
     <OTHER ARGUMENTS>) {}

2. Modify the internal C++ function in src/model_*.cpp to accept the mortality rate as an argument of the type const double; this also applies to integer-like inputs as the double type is better suited for R’s numeric type, which most users will be working with. This is level 2 of the code, one level down.

   // Note this function is exposed to R but not exported from the package
   //' [[Rcpp::export(name=".model_NAME_cpp")]]    // <== Use model name.
   Rcpp::List model_NAME_internal(
     const Eigen::MatrixXd &initial_state,
     const double &transmission_rate,
     const double &infectiousness_rate,
     const double &recovery_rate,
     const double &mortality_rate,       // <== New parameter added.
     <OTHER_ARGUMENTS>) {
       /* model setup and ODE solving */
     }

3. Add the mortality rate to the std::unordered_map of key-value pairs of parameters and their names, in src/model_*.cpp; this is used to access the parameters in the ODE solver, as well as enabling rate interventions and time-dependence.

   Rcpp::List model_NAME_internal(
     const Eigen::MatrixXd &initial_state,
     const double &transmission_rate,
     const double &infectiousness_rate,
     const double &recovery_rate,
     const double &mortality_rate
     <OTHER_ARGUMENTS>) {
   
       // create a map of the model parameters
       std::unordered_map<std::string, double> model_params{
           {"transmission_rate", transmission_rate},
           {"infectiousness_rate", infectiousness_rate},
           {"recovery_rate", recovery_rate},
           {"mortality_rate", mortality_rate}    // <== New parameter added.
         };
   
       /* ODE solving */
     }

4. Modify the compartmental transitions in the FunctionObject ODE model inst/include/model_*.h appropriately to include the mortality rate. This is level 3 of the code, two levels down.

   namespace epidemics {
   
   struct epidemic_new {
     // Model elements as struct members
     // two maps for the model parameters, one for dynamic modification
     const std::unordered_map<std::string, double> model_params;
     std::unordered_map<std::string, double> model_params_temp;
     <OTHER STRUCT MEMBERS>
   
     // constructor
     epidemic_new(<ARGUMENTS>) : <CONSTRUCTOR LIST> {}
   
     // overloaded operator
     void operator()(const odetools::state_type& x,    // the current state
               odetools::state_type& dxdt,             // change in state
               const double t) {                       // the current time
       /* dynamic parameter calculation here */
   
       // implement background mortality rate here
       // NOTE: compartments are hard-coded as column positions
       // and demographic groups are rows.
       // Accessed columns must be converted to arrays for
       // vectorised operations.
   
       // For the 0-th column, i.e., susceptibles
       dxdt.col(0) = <NEW INFECTIONS AND VACCINATIONS> -
         (
           model_params_temp.at("mortality_rate") *    // mortality rate *
             x.col(0).array()                          // susceptibles
         );
   
       // and so on for all i columns in x
     }
   
   }

5. Write or update unit tests to check for the statistical correctness of the model with the newly added parameter (e.g., for this example, test that a non-zero mortality rate leads to fewer individuals at the end of the model than at the start).

Modifying compartmental flows without changing compartments

Modifying compartmental flows is straightforward for the ODE models in epidemics.

This section focuses on adding novel inflows and outflows via some illustrative examples. The initial steps are similar to those for adding a model parameter, and are summarised here. The section on adding group-specific parameters covers how to add parameters with dimensions > 1.

1. First, navigate to the end of each model FunctionObject defined in inst/include/model_*.h and edit the ODEs found there (see examples below). Users can hard-code the parameter values in these files for quick edits, but we recommend taking the next few steps to allow changing these values from R.

2. Next, modify the C++ function arguments in src/model_*.cpp so that any newly added rate parameters can be passed from R to C++. Typically, single parameters should be passed by reference as const double.

• If you seek to use time-dependence and rate interventions with the newly added parameters should follow the steps in the section above.

3. Finally, modify the R code in R/model_*.R so that the R function accepts the new parameter as an appropriately checked argument.

4. When adding parameters such as birth and death rates, make sure that this is accounted for in any unit tests (see especially the current expectation in most models that there is a constant population size).

// define a uniform spillover rate from contact with animals
const double spillover_rate = SOME_VALUE;

// include this additional S -> E transition
// note x.col(0).array() is the vector of age-specific susceptibles
// each model header file 
const Eigen::ArrayXd zoonotic_infection = spillover_rate * x.col(0).array();

dxdt.col(0) = -sToE - sToV - zoonotic_infection;  // -β*S*contacts*I - ν - zoonotic infection

// define an age-specific waning rate
const double waning_rate = SOME_VALUE;

// include this new R -> S transition
// note x.col(3).array() is the vector of age-specific recovereds
const Eigen::ArrayXd re_susceptibles = waning_rate * x.col(3).array();

// add individuals with renewed susceptibility
dxdt.col(0) = -sToE - sToV + re_susceptibles;  // -β*S*contacts*I - ν - re_susceptibles

Modifying which parameters can be time dependent

Making the infection parameters of existing models time-dependent is straightforward using the time_dependence argument to model functions. Currently, all models support time-dependence on most infection parameters. This does not include the vaccination rate (see more on that below).

See the vignette on time-dependence for how this functionality currently works.

To modify which parameters of new or existing models can be made (optionally) time-dependent using the existing functionality:

1. For ODE models, first modify the std::unordered_map named model_params in the corresponding model source file, src/model_*.cpp, to add or remove a key-value pair. Adding a parameter to this map automatically enables both time-dependence as well as rate interventions; removing a parameter from this map means disallowing both.

2. Next, modify the model functions and the argument cross-checking functions to ensure that the model allows time-dependence and rate interventions on the correct set of infection parameters.

   • The vector of allowed targets for time-dependence is an argument to the function .cross_check_timedep() called in model_*() from R/model_*.R, while the vector of allowed targets for rate interventions is a similar argument to the function .cross_check_intervention() called in .check_prepare_args_*() from R/check_args_*.R; in all cases * refers to the model name.

   • This initial cross-checking is required to prevent errors in the C++ code (typically, a key-not-found error) which may be difficult for users to understand and fix.

3. Finally, ODE model code can be simplified to remove time-dependence (and rate interventions, or one of them) entirely by deleting the call to time_dependence::apply_time_dependence in the ODE model FunctionObject in inst/include/model_*.h. This will not affect upstream code, and may be one way of gaining a small improvement in model performance, so long as time-dependence is not required.

4. For the stochastic Ebola model, infection parameters are held in a named list called parameters. The elements of this list that are targeted for time-dependence are modified within the model loop.

   • Adding parameters to the list allows them to be be targeted by inbuilt time-dependence functionality.

   • However, removing parameters from the list will lead to errors as they are needed for compartmental transitions. Instead, add input checks on the arguments intervention and time_dependence to prevent these elements from being modified.

   • Remove the appropriate code from within the loop to completely disallow time-dependence and rate interventions without affecting upstream code.

   • Note that the Ebola model does support modifying the infectiousness rate or the removal rate, as these are used to calculate the number of sub-compartments (representing days) in each epidemiological compartment; these values must remain fixed throughout the model’s run time.

Group-specific infection parameters

epidemics ODE models can be modified to support group-specific infection parameters instead of uniform, population-wide parameters. This is essentially how vaccinations are implemented. Note that the Ebola model does not support age-stratification and is not discussed here.

Since epidemics focuses on general models for use within the first 100 days of a pandemic, it does not currently support groups-specific parameters. Supporting these parameters in the ODE models would require a substantial rewrite of how epidemics is vectorised.

1. First, in the C++ function signature in src/model_*.cpp change the type of the parameter to Eigen::ArrayXd.

2. Next, in src/model_*.cpp, remove the parameter from inclusion in the std::unordered_map of parameters to prevent initialisation errors (as the map types are std::string and double). This means time-dependence and rate interventions cannot be applied to this parameter.

3. Modify the model ODE struct constructor in inst/include/model_*.h to require an Eigen::ArrayXd for the group-specific parameter, and define the parameter as a struct member (initialised from the constructor).

4. Modify ODEs to use the parameter vector/array as appropriate, by simply replacing the single parameter with the parameter vector in the ODE code.

5. Ensure that all R wrapper code – especially input checking and cross-checking code – is updated to reflect that one or more infection parameters are now vectors.

Adding epidemiological compartments

Adding epidemiological compartments is similar to modifying compartmental transitions, but requires a few extra steps.

1. Add and modify model ODEs as appropriate in inst/include/model_*.h (see sections above).

2. Modify the compartments vector in R/model_*.R to include the name of the new compartment. This is used both in checking that the input <population> has the correct number of compartments in its initial state matrix, as well as for formatting the output data.

Changing vaccination rates over time

epidemics treats vaccination as special process that requires a composable element in the form of an S3 class, <vaccination>, which allows defining time-limited, group-specific vaccination rates for multiple vaccination doses if necessary.

Vaccination is only strictly allowed in the ODE models ‘default’ and ‘Vacamole’, and these are covered here. Begin by dismantling the model code to remove any functionality related to the <vaccination> class.

1. Modify R/model_*.R to remove vaccination from the function arguments, and any input checking on the vaccination argument. Also remove the inclusion of vaccination in the model arguments <data.table>.

2. Remove any cross-checking of the vaccination against the population in R/check_args_*.R, and remove the inclusion of vax_time_begin, vax_time_end, and vax_nu from the return of the function .check_prepare_args_*().

3. Remove the vax_time_begin, vax_time_end, and vax_nu arguments from the model ODE struct constructor in inst/include/model_*.h; also remove these from the model object initialisation in src/model_*.cpp. Remove any code related to calculating current vaccination rates from the model ODE code.

4. At this stage, the ODE model does not support vaccination using a <vaccination>, and any existing “vaccinated” compartment is essentially isolated from the other compartments, as the value of current_nu

5. Next, follow the steps described above to re-establish compartmental transitions from the appropriate compartments to the vaccinated compartment, using either a single population-wide vaccination rate, or an age-specific rate.

6. Note that the number of individuals vaccinated should be calculated as counts, and not a rate, when used in compartmental transitions. This is because the rate of vaccination is typically specified as a proportion of the total population, or of age groups, and should be dependent on the number of individuals available to vaccinate.