Default model with odin-backend example

[pratikunterwegs]

Followin from a meeting with @adamkucharski @sbfnk and @MJomaba and @LloydChapman, this issue is to request an implementation of the default model using {odin} as a backend. Having the implementation as a vignette initially could be a good starting point and help explore how the two packages could be integrted.

We'd want to see the current epidemics features preserved:

1. Age stratification with age-dependent social contacts;
2. Time-dependence on infection parameters;
3. Interventions on social contacts, and vaccination specified by a <vaccination>.

[pratikunterwegs]

Here's a brief initial reprex for an SIR model with age-stratified social contacts, and a crude intervention on the transmission rate. Proof of concept only, but SEIR is only a few changes away.

Some thoughts:

1. Usability could be improved by having some way of passing the initial conditions more flexibly, as they are (have to be?) hardcoded currently - will look into this;
2. I think it should be possible to pass time-limitation flags as 2D or 3D arrays to implement interventions and vaccination regimes, but passing the {epidemics} classes themselves might be iffy.

library(odin)
library(socialmixr)
#> 
#> Attaching package: 'socialmixr'
#> The following object is masked from 'package:utils':
#> 
#>     cite
library(data.table)
library(ggplot2)

# get social contacts data
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 20, 40),
  symmetric = TRUE
)
#> Using POLYMOD social contact data. To cite this in a publication, use the 'cite' function
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option

# assume arbitrary populations for each demo group
population_size <- rep(1e6, 3)

# prepare contact matrix, divide by leading eigenvalue and rowwise by popsize
contact_matrix <- t(contact_data[["matrix"]])
contact_matrix <- (contact_matrix / max(Re(eigen(contact_matrix)$values))) /
  population_size

# View
contact_matrix
#>                  
#> contact.age.group         [,1]         [,2]         [,3]
#>           [0,20)  6.681250e-07 2.367991e-07 1.326870e-07
#>           [20,40) 2.644325e-07 4.114381e-07 2.224524e-07
#>           40+     2.596591e-07 3.898319e-07 4.242123e-07

# define SIR model
sir <- odin::odin({
  # NOTE: crude time-dependence of transmission rate
  # similar to a `<rate_intervention>`
  beta_ <- if (t > 210 && t < 300) beta / 2 else beta

  # number of age groups taken from contact matrix passed by user
  n <- user()

  # FOI is contacts * infectious * transmission rate
  # returns a matrix, must be converted to a vector/1D array
  lambda_prod[, ] <- C[i, j] * I[j] * beta_
  lambda[] <- sum(lambda_prod[i, ])

  ## Derivatives
  deriv(S[]) <- -S[i] * lambda[i]
  deriv(I[]) <- (S[i] * lambda[i]) - (gamma * I[i])
  deriv(R[]) <- gamma * I[i]

  ## Initial conditions: seems these cannot be passed by user?
  initial(S[]) <- 1e6 - 1
  initial(I[]) <- 1
  initial(R[]) <- 0

  ## parameters
  beta <- 1.3 / 7
  gamma <- 1 / 7
  C[, ] <- user() # user defined contact matrix

  # dimensions - all rely on contact matrix
  dim(lambda_prod) <- c(n, n)
  dim(lambda) <- n
  dim(S) <- n
  dim(I) <- n
  dim(R) <- n
  dim(C) <- c(n, n)
})
#> Loading required namespace: pkgbuild
#> Generating model in c
#> ℹ Re-compiling odin955a4062 (debug build)
#> ── R CMD INSTALL ───────────────────────────────────────────────────────────────
#> * installing *source* package ‘odin955a4062’ ...
#> ** using staged installation
#> ** libs
#> using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
#> using SDK: ‘MacOSX14.0.sdk’
#> clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -UNDEBUG -Wall -pedantic -g -O0 -c odin.c -o odin.o
#> clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -UNDEBUG -Wall -pedantic -g -O0 -c registration.c -o registration.o
#> clang -arch arm64 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/Library/Frameworks/R.framework/Resources/lib -L/opt/R/arm64/lib -o odin955a4062.so odin.o registration.o -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation
#> ld: warning: -single_module is obsolete
#> ld: warning: -multiply_defined is obsolete
#> installing to /private/var/folders/3t/nk820_1d4fvf8r90j7j1ldbr0000gp/T/RtmpKSY1jr/devtools_install_70bc1415b97d/00LOCK-file70bc4d32de76/00new/odin955a4062/libs
#> ** checking absolute paths in shared objects and dynamic libraries
#> * DONE (odin955a4062)
#> ℹ Loading odin955a4062

# initialise model and run over time 0 - 600
mod <- sir$new(C = contact_matrix, n = nrow(contact_matrix))
t <- seq(0, 600)
y <- mod$run(t)

# convert to data.table and plot infectious in each age class
y <- as.data.table(y)
y <- melt(y, id.vars = c("t"))

ggplot(y[variable %like% "I"]) +
  geom_line(aes(t, value, col = variable)) +
  facet_grid(~variable)

^{Created on 2024-04-25 with reprex v2.0.2}

[adamkucharski]

Thanks for laying this out. It should be possible to pass input parameters, e.g. see early code in this vignette @LloydChapman put together and corresponding odin model.

Perhaps a starting point would be to implement the 'getting started' {epidemics} model as {odin} (or a similar vignette if that's easier/quicker), so can check outputs all match, then identify places we could streamline or improve integration?

[pratikunterwegs]

Thanks, got it. One thing that's seem to be more difficult than I initially thought is event handling of time-limited interventions. As in the example, it's relatively easy to pass a single uniform intervention within a single range, but I'm unsure whether the DSL allows for multiple overlapping interventions affecting age groups differently (which we'd want for policy-relevant scenarios e.g. school/work closures). Happy to get any pointers on this from @LloydChapman!

[adamkucharski]

Thanks. So sounds like it should be feasible to reproduce the basic README example while we dig more into how to implement overlapping interventions?

It must be doable in theory given {sircovid} etc. (I'm not sure how much overlap we had with interventions in FP analysis). But, of course, depends how much of this needs to be hard coded in DSL, versus the flexibility of defining as inputs for current {epidemics} models

[pratikunterwegs]

Yes, that's the example shown as a reprex here already. Since @avallecam is also referencing this issue might it be better to add this to the How To repo instead? Would also save on adding dependencies.

If I'm right sircovid implemented interventions as time-dependence on transmission rate, see linked issue and related project(as in the reprex) rather than group specific contact reduction. That might be doable but equally I think it's where epidemics steals a march.

[adamkucharski]

Yep, howto seems like good home for now. For full consistency, would be useful to see if it's possible for above reprex to be wrapped into a function call for existing script (rather than needing ICs, parameters etc. to be defined in the odin code), e.g. so below function calls a compiled odin backend, rather than current one:

model_default_odin(
  population = uk_population,
  transmission_rate = 1.5 / 7.0,
  infectiousness_rate = 1.0 / 3.0,
  recovery_rate = 1.0 / 7.0,
  intervention = list(contacts = close_schools),
  time_end = 600, increment = 1.0
)

[pratikunterwegs]

Sounds good, I'll figure out the wrapping and make a PR to {howto}.

[LloydChapman]

Thanks, got it. One thing that's seem to be more difficult than I initially thought is event handling of time-limited interventions. As in the example, it's relatively easy to pass a single uniform intervention within a single range, but I'm unsure whether the DSL allows for multiple overlapping interventions affecting age groups differently (which we'd want for policy-relevant scenarios e.g. school/work closures). Happy to get any pointers on this from @LloydChapman!

Good question. From looking through the Imperial papers, it looks like they never actually modelled interventions like school closures as affecting different age groups differently (they just changed the overall transmission rate coefficient). It should be perfectly possible to do, though, by making the "contact" matrix (or strictly speaking "transmission" matrix) in the code a time-dependent function of some user inputs, but then you do need to hard code this in the model. There's some flexibility in terms of how many changepoints you have and how "parameter" values change, as you can specify a vector/list of input values for the parameter/contact matrix at the changepoints and then convert these to a time-dependent function as e.g. here and here, where I'm pretty sure the function that does the conversion can be anything that interpolates the values at the changepoints over all the simulated time points.

[pratikunterwegs]

Thanks @LloydChapman, I'll take a look but as I'm not an {odin} user the linked examples might be difficult to understand right away - could you provide a reprex of what you mean? Just to clarify, does {odin} support matrix operations and function definitions within the model code?

[sbfnk]

  ## Initial conditions: seems these cannot be passed by user?

I think you can just define e.g. a user-defined matrix and then use that, i.e.

  init[, ] <- user()
  dim(init) <- c(3, n)
  initial(S[]) <- init[1, i]
  initial(I[]) <- init[2, i]
  initial(R[]) <- init[3, i]

and then in the call

mod <- sir$new(
  C = contact_matrix,
  n = nrow(contact_matrix),
  init = matrix(c(rep(1e6 - 1, 3), rep(1, 3), rep(0, 3)), byrow = TRUE, ncol = 3)
)

[sbfnk]

Thanks @LloydChapman, I'll take a look but as I'm not an {odin} user the linked examples might be difficult to understand right away - could you provide a reprex of what you mean? Just to clarify, does {odin} support matrix operations and function definitions within the model code?

Don't know much about odin but looking at your existing code could you do something like the following?

  C[, ] <- user() # user defined contact matrix
  C_[, ] <- C[i, j]
  C_[1, 1] <- if (t > 210 && t < 300) C[1, 1] / 2 else C[1, 1]

[adamkucharski]

The {seir} package vignette has also been updated (thanks @LloydChapman for also fixing the fitting issue!) so hopefully clearer here how to pass pre-defined inputs into model: https://github.com/LloydChapman/seir/blob/main/vignettes/fit_seirdage.Rmd