R for modellers - Vignette 11

Solving DDE in R

Julien Arino

Department of Mathematics

University of Manitoba*

* The University of Manitoba campuses are located on original lands of Anishinaabeg, Cree, Oji-Cree, Dakota and Dene peoples, and on the homeland of the Métis Nation.

Outline

• The deSolve library
• Example - Fitting data

The deSolve library

• As I have already pointed out, deSolve: > The functions provide an interface to the FORTRAN functions ‘lsoda’, ‘lsodar’, ‘lsode’, ‘lsodes’ of the ‘ODEPACK’ collection, to the FORTRAN functions ‘dvode’, ‘zvode’ and ‘daspk’ and a C-implementation of solvers of the ‘Runge-Kutta’ family with fixed or variable time steps

• So you are benefiting from years and year of experience: ODEPACK is a set of Fortran (77!) solvers developed at Lawrence Livermore National Laboratory (LLNL) starting in the late 70s

• Other good solvers are also included, those written in C

• Refer to the package help for details

Using deSolve for simple ODEs

As with more numerical solvers, you need to write a function returning the value of the right hand side of your equation (the vector field) at a given point in phase space, then call this function from the solver

library(deSolve)
rhs_logistic <- function(t, x, p) {
  with(as.list(x), {
    dN <- p$r * N *(1-N/p$K)
    return(list(dN))
  })
}
params = list(r = 0.1, K = 100)
IC = c(N = 50)
times = seq(0, 100, 1)
sol <- ode(IC, times, rhs_logistic, params)

This also works: add p to arguments of as.list and thus use without p$ prefix

library(deSolve)
rhs_logistic <- function(t, x, p) {
  with(as.list(c(x, p)), {
    dN <- r * N *(1-N/K)
    return(list(dN))
  })
}
params = list(r = 0.1, K = 100)
IC = c(N = 50)
times = seq(0, 100, 1)
sol <- ode(IC, times, rhs_logistic, params)

In this case, beware of not having a variable and a parameter with the same name..

Default method: lsoda

• lsoda switches automatically between stiff and nonstiff methods

• You can also specify other methods: “lsode”, “lsodes”, “lsodar”, “vode”, “daspk”, “euler”, “rk4”, “ode23”, “ode45”, “radau”, “bdf”, “bdf_d”, “adams”, “impAdams” or “impAdams_d” ,“iteration” (the latter for discrete-time systems)

ode(y, times, func, parms, 
    method = "ode45")

• You can even implement your own integration method

Example - Fitting data

Example - Fitting data

• Note that this is a super simplified version of what to do
• Much more elaborate procedures exist
  • Roda. Bayesian inference for dynamical systems
  • Portet. A primer on model selection using the Akaike Information Criterion
• Let us grab some epi data online and fit an SIR model to it
• Don’t expect anything funky, as I said, this is the baby version
• Also, keep in mind that any identification procedure is subject to risks due to identifiability issues; see, e.g., Roda et al, Why is it difficult to accurately predict the COVID-19 epidemic?

Principle

• Data is a set \((t_i,y_i)\), \(i=1,\ldots,N\), where \(t_i\in\mathcal{I}\), some interval
• Solution to SIR is \((t,x(t))\) for \(t\in\mathcal{I}\)
• Suppose parameters of the model are \(p\)
• We want to minimise the error function \[ E(p) = \sum_{i=1}^N \|x(t_i)-y_i\| \]
• Norm is typically Euclidean, but could be different depending on objectives
• So given a point \(p\) in (admissible) parameter space, we compute the solution to the ODE, compute \(E(p)\)
• Using some minimisation algorithm, we seek a minimum of \(E(p)\) by varying \(p\)

What are \(y_i\) and \(x(t_i)\) here?

• In epi data for infectious diseases, we typically have incidence, i.e., number of new cases per unit time
• In SIR model, this is \(\beta SI\) or \(\beta SI/N\), so, if using mass action incidence and Euclidean norm \[ E(p)=\sum_{i=1}^N(\beta S(t_i)I(t_i)-y_i)^2 \] or, if using standard incidence \[ E(p)=\sum_{i=1}^N \left(\beta \frac{S(t_i)I(t_i)}{N}-y_i\right)^2 \]

Implementing in practice

See the code practicum_01_fitting.R, which we will go over now