library(epidemics)
library(socialmixr)
library(odin)
library(data.table)
library(ggplot2)How to use epidemics classes with odin-generated models
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.
Steps in code
Prepare population initial conditions
# get social contacts data
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
polymod,
countries = "United Kingdom",
age.limits = c(0, 20, 40),
symmetric = TRUE
)
#> 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
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
#> <contacts_intervention> object
#>
#> Intervention name:
#> NA
#>
#> Begins at:
#> [1] 200
#>
#> Ends at:
#> [1] 260
#>
#> Reduction:
#> Interv. 1
#> Demo. grp. 1 0.50
#> Demo. grp. 2 0.01
#> Demo. grp. 3 0.01
# 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_vectorDefine an epidemic model in odin
This is an SIR model that omits the ‘exposed’ compartment for this demonstration.
# 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
})
#> Loading required namespace: pkgbuild
#> Generating model in c
#> ℹ Re-compiling odin6501fa20 (debug build)
#> ── R CMD INSTALL ───────────────────────────────────────────────────────────────
#> * installing *source* package ‘odin6501fa20’ ...
#> ** using staged installation
#> ** libs
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> gcc -I"/opt/R/4.4.2/lib/R/include" -DNDEBUG -I/usr/local/include -fpic -g -O2 -UNDEBUG -Wall -pedantic -g -O0 -c odin.c -o odin.o
#> odin.c: In function ‘user_get_scalar_int’:
#> odin.c:398:47: warning: format ‘%d’ expects argument of type ‘int’, but argument 2 has type ‘const char *’ [-Wformat=]
#> 398 | Rf_error("Expected scalar integer for '%d'", name);
#> | ~^ ~~~~
#> | | |
#> | int const char *
#> | %s
#> odin.c: In function ‘user_get_array’:
#> odin.c:556:48: warning: format ‘%d’ expects argument of type ‘int’, but argument 2 has type ‘size_t’ {aka ‘long unsigned int’} [-Wformat=]
#> 556 | Rf_error("Incorrect size of dimension %d of %s (expected %d)",
#> | ~^
#> | |
#> | int
#> | %ld
#> 557 | i + 1, name, dim_expected);
#> | ~~~~~
#> | |
#> | size_t {aka long unsigned int}
#> gcc -I"/opt/R/4.4.2/lib/R/include" -DNDEBUG -I/usr/local/include -fpic -g -O2 -UNDEBUG -Wall -pedantic -g -O0 -c registration.c -o registration.o
#> gcc -shared -L/opt/R/4.4.2/lib/R/lib -L/usr/local/lib -o odin6501fa20.so odin.o registration.o -L/opt/R/4.4.2/lib/R/lib -lR
#> installing to /tmp/Rtmpih3Mpy/devtools_install_17e23602140f/00LOCK-file17e24b519637/00new/odin6501fa20/libs
#> ** checking absolute paths in shared objects and dynamic libraries
#> * DONE (odin6501fa20)
#> ℹ Loading odin6501fa20# 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)