Workshop 1: calibrating a deterministic model

Danny Scarponi, Andy Iskauskas

1 Objectives

This workshop offers an interactive introduction to the main functionality of the hmer package, using a simple example of an epidemiological model. hmer allows you to efficiently implement the Bayes Linear emulation and history matching process, a calibration method that has been successfully employed in several sciences (epidemiology, cosmology, climate science, systems biology, geology, energy systems etc.). Note that, even though this workshop focuses on an epidemiological model, hmer can be used to calibrate complex models arising in any field.

In this workshop, you will be invited to carry out a series of tasks (see “Task” boxes) which will enhance your understanding of the package and its tools. Thanks to these activities, you will learn to calibrate deterministic models using history matching and model emulation, and to use your judgement to customise the process. This workshop should be considered as a natural continuation of Tutorial 1, which introduces the history matching with emulation framework with a one-dimensional example, and of Tutorial 2, which gives a general overview of the history matching with emulation process for deterministic models, and shows how to perform it using hmer. Following Workshop 1, you may also want to read Workshop 2, where we demonstrate how to calibrate stochastic models.

For further discussion and justification of the various stages of the history matching with emulation process, please see Bower, Goldstein, and Vernon (2010) and Vernon et al. (2018).

Note that when running the workshop code on your device, you should not expect the hmer visualisation tools to produce the same exact output as the one you can find in the following sections. This is mainly because the maximinLHS function, that you will use to define the initial parameter sets on which emulators are trained, does return different Latin Hypercube designs at each call.

Before starting the tutorial, you will need to run the code contained in the box below: it will load all relevant dependencies and a few helper functions which will be introduced and used later.

Show: Code to load relevant libraries and helper functions

library(hmer)
library(deSolve)
library(ggplot2)
library(reshape2)
library(purrr)
library(tidyverse)
library(lhs)
set.seed(123)
############################# HELPER FUNCTIONS #############################
# `ode_results` provides us with the solution of the differential equations for a given 
# set of parameters. This function assumes an initial population of 
# 900 susceptible individuals, 100 exposed individuals, and no infectious 
# or recovered individuals.
ode_results <- function(parms, end_time = 365*2) {
  forcer = matrix(c(0, parms['beta1'], 100, parms['beta2'], 180, parms['beta2'], 270, parms['beta3']), 
                  ncol = 2, byrow = TRUE)
  force_func = approxfun(x = forcer[,1], y = forcer[,2], method = "linear", rule = 2)
  des = function(time, state, parms) {
    with(as.list(c(state, parms)), {
      dS <- b*(S+E+I+R)-force_func(time)*I*S/(S+E+I+R) + omega*R - mu*S
      dE <- force_func(time)*I*S/(S+E+I+R) - epsilon*E - mu*E
      dI <- epsilon*E - alpha*I - gamma*I - mu*I
      dR <- gamma*I - omega*R - mu*R
      return(list(c(dS, dE, dI, dR)))
    })
  }
  yini = c(S = 900, E = 100, I = 0, R = 0)
  times = seq(0, end_time, by = 1)
  out = deSolve::ode(yini, times, des, parms)
  return(out)
}
# `get_results` acts as `ode_results`, but has the additional feature 
# of allowing us to decide which outputs and times should be returned. 
# For example, to obtain the number of infected and susceptible individuals 
# at t=25 and t=50, we would set `times=c(25,50)` and `outputs=c('I','S')`.
get_results <- function(params, times, outputs) {
  t_max <- max(times)
  all_res <- ode_results(params, t_max)
  actual_res <- all_res[all_res[,'time'] %in% times, c('time', outputs)]
  shaped <- reshape2::melt(actual_res[,outputs])
  return(setNames(shaped$value, paste0(shaped$Var2, actual_res[,'time'], sep = "")))
}

References

Bower, Richard G, Michael Goldstein, and Ian Vernon. 2010. “Galaxy Formation: A Bayesian Uncertainty Analysis.” Bayesian Analysis 5 (4): 619–69.

Vernon, Ian, Junli Liu, Michael Goldstein, James Rowe, Jen Topping, and Keith Lindsey. 2018. “Bayesian Uncertainty Analysis for Complex Systems Biology Models: Emulation, Global Parameter Searches and Evaluation of Gene Functions.” BMC Systems Biology 12 (1): 1–29.