How to use epidemics-classes with {odin} models

analyses/simulate_transmission/epidemics_odin_interop.qmd

@@ -0,0 +1,153 @@
+---
+title: "How to use _epidemics_ classes with _odin_-generated models"
+format: 
+  html:
+    code-link: true
+editor: source
+editor_options: 
+  chunk_output_type: console
+date: last-modified
+toc: true
+toc_float: true
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+## Ingredients
+
+This guide shows how to use some of the convenient features of _epidemics_, especially the classes that describe populations, interventions, and vaccination regimes, with models that are generated by _odin_.
+
+This guide shows this interoperability by implementing the default model from _epidemics_ in _odin_.
+
+```{r}
+library(epidemics)
+library(socialmixr)
+library(odin)
+library(data.table)
+library(ggplot2)
+```
+
+## Steps in code
+
+### Prepare population initial conditions
+
+```{r}
+# get social contacts data
+polymod <- socialmixr::polymod
+contact_data <- socialmixr::contact_matrix(
+  polymod,
+  countries = "United Kingdom",
+  age.limits = c(0, 20, 40),
+  symmetric = TRUE
+)
+
+# assume arbitrary populations for each demo group
+demography_vector <- contact_data$demography$population
+
+# 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))) /
+  demography_vector
+
+# an intervention to close schools
+close_schools <- intervention(
+  type = "contacts",
+  time_begin = 200,
+  time_end = 260,
+  reduction = matrix(c(0.5, 0.01, 0.01))
+)
+
+# view the intervention
+close_schools
+
+# initial conditions defined as proportions and converted to absolute values
+initial_conditions <- matrix(
+  c(S = 1 - 1e-6, I = 1e-6, R = 0),
+  nrow = length(demography_vector),
+  ncol = 3L, byrow = TRUE
+) * demography_vector
+```
+
+### Define an epidemic model in _odin_
+
+This is an SIR model that omits the 'exposed' compartment for this demonstration.
+
+```{r}
+# define SIR model
+sir <- odin::odin({
+  # NOTE: crude time-dependence of transmission rate
+  # similar to a `<rate_intervention>`
+  beta_t[] <- if (t > intervention_interval[1] && t < intervention_interval[2]) beta * (1.0 - reduction[i]) else beta*(1 + 0*reduction[i])
+
+  # 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_t[j]
+  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: passed by user
+  initial(S[]) <- init_S[i]
+  initial(I[]) <- init_I[i]
+  initial(R[]) <- init_R[i]
+
+  ## parameters
+  beta <- user()
+  gamma <- 1 / 7
+  C[, ] <- user() # user defined contact matrix
+  init_S[] <- user()
+  init_I[] <- user()
+  init_R[] <- user()
+  reduction[] <- user()
+  intervention_interval[] <- user()
+
+  # 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(init_S) <- n
+  dim(init_I) <- n
+  dim(init_R) <- n
+  dim(reduction) <- n
+  dim(beta_t) <- n
+  dim(C) <- c(n, n)
+  dim(intervention_interval) <- 2
+})
+```
+
+```{r}
+# initialise model and run over time 0 - 600
+mod <- sir$new(
+  beta = 1.3 / 7,
+  reduction = close_schools$reduction[,],
+  intervention_interval = c(close_schools$time_begin, close_schools$time_end),
+  C = contact_matrix, n = nrow(contact_matrix),
+  init_S = initial_conditions[, 1],
+  init_I = initial_conditions[, 2],
+  init_R = initial_conditions[, 3]
+)
+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)) +
+  labs(x = "Time", y = "# Infectious", col = "Age group") +
+  facet_grid(~variable)
+```