Additional Learning

Here we provide examples of epidemiological models, visualizations, and simulation strategies using epiworldR. These examples are intended to serve as jumping off points for new material and as quick references to the material covered in the workshop, thus we often simply show the code with little to no explanation. For further learning, see the workshop Parts 1-3 or the epiworldR package documentation.

Epidemiological Model Examples

Examples of popular epidemiological models implemented in epiworldR.

SIR Model (Fully Connected Graph)

Code
model_sirconn <- ModelSIRCONN(
  name                = "COVID-19",
  n                   = 10000,
  prevalence          = 0.01,
  contact_rate        = 5,
  transmission_rate   = 0.4,
  recovery_rate       = 0.95
  ) |>
  verbose_off() |>
  run(ndays = 50, seed = 1912) |>
  plot()

SIR Model (Random Graph)

Code
model_sir <- ModelSIR(
  name              = "COVID-19",
  prevalence        = .01,
  transmission_rate = .7,
  recovery_rate     =  .3
  ) |>
  agents_smallworld(n = 100000, k = 10, d = FALSE, p = .01) |>
  verbose_off() |>
  run(ndays = 50, seed = 1912) |>
  plot()

SEIR Model (Fully Connected Graph)

Code
model_seirconn <- ModelSEIRCONN(
  name                = "COVID-19",
  prevalence          = 0.01, 
  n                   = 10000,
  contact_rate        = 4, 
  incubation_days     = 7, 
  transmission_rate   = 0.6,
  recovery_rate       = 0.5
  ) |>
  verbose_off() |>
  run(ndays = 50, seed = 1912) |>
  plot()

SEIR Model (Random Graph)

Code
model_seir <- ModelSEIR(
  name              = "COVID-19", 
  prevalence        = 0.01,
  transmission_rate = 0.9, 
  recovery_rate     = 0.1,
  incubation_days   = 4
  ) |>
  agents_smallworld(n = 1000, k = 5, d = FALSE, p = .01) |>
  verbose_off() |>
  run(ndays = 50, seed = 1912) |>
  plot()

SIR Logit

Code
set.seed(2223)
n <- 100000

X <- cbind(
  Intercept = 1,
  Female    = sample.int(2, n, replace = TRUE) - 1
  )

coef_infect  <- c(.1, -2, 2)
coef_recover <- rnorm(2)

model_logit <- ModelSIRLogit(
  vname             = "COVID-19",
  data              = X,
  coefs_infect      = coef_infect,
  coefs_recover     = coef_recover, 
  coef_infect_cols  = 1L:ncol(X),
  coef_recover_cols = 1L:ncol(X), 
  prob_infection    = .8,
  recovery_rate     = .3,
  prevalence        = .01
  ) |>
  agents_smallworld(n = n, k = 8, d = FALSE, p = .01) |>
  verbose_off() |>
  run(ndays = 50, seed = 1912) |>
  plot()

Example Scenario: Comorbidities Using Logit Functions

Often, we want to model the effects of comorbidities on a disease. In this example, we’ll examine the effects of obesity on the probability of recovery from the flu.

Model Setup

Create two identical models using the ModelSEIRCONN() function. One will have comorbidities, the other will not.

Code
# With comorbidities
model_comorbid <- ModelSEIRCONN(
  name              = "Flu",
  n                 = 10000, 
  prevalence        = 0.001, 
  contact_rate      = 2.1,
  transmission_rate = 0.5,
  incubation_days   = 7,
  recovery_rate     = 1/4
  )

# Without comorbidities
model_no_comorbid <- ModelSEIRCONN(
  name              = "Flu",
  n                 = 10000, 
  prevalence        = 0.001, 
  contact_rate      = 2.1,
  transmission_rate = 0.5,
  incubation_days   = 7,
  recovery_rate     = 1/4
  )

Add Comorbidities

The file part2b_comorb.rds contains obesity data for the agents. This is formatted as a matrix with two columns: baseline and obesity. We’ll read this in and assign it to the agents of the comorbidity model using the set_agents_data() function.

Code
obesity_data <- readRDS("part2b_comorb.rds")
set_agents_data(model_comorbid, obesity_data)

Now that population includes some agents with the obesity condition.

Set Recovery Rate Based on Comorbidity

Use the virus_fun_logit() function to create a function for the probability of recovery. The function takes the following parameters:

  • vars: the variables to use in the model
  • coefs: the coefficients for each variable
  • model: the model object
Code
# Logit function 
lfun <- virus_fun_logit(
  vars  = 0:1,
  coefs = c(-1.0986, -0.8472), 
  model_comorbid
  )
Note

To build the logit function, we used the following: (a) Under the logit model, the coefficient needed for the baseline probability of 0.25 is computed using qlogis(0.25). With that, we can compute the associated coefficient to obese individuals with plogis(qlogis(.25) + x) = .125 -> qlogis(.25) + x = plogis(.125) -> x = qlogis(.125) - qlogis(.25)

Use the set_prob_recovery_fun() function to set the probability of recovery function for the virus to the logit function.

Code
set_prob_recovery_fun(
  virus = get_virus(model_comorbid, 0), 
  model = model_comorbid,
  vfun  = lfun
  )

Run the Model

Run both models for the 50 days with the same random seed and compare the results.

Code
verbose_off(model_comorbid)
verbose_off(model_no_comorbid)
run(model_comorbid, ndays = 50, seed = 1231)
run(model_no_comorbid, ndays = 50, seed = 1231)

op <- par(mfrow = c(1, 2), cex = .7)
plot_incidence(model_no_comorbid, main = "Without Comorbidities")
plot_incidence(model_comorbid, main = "With Comorbidities")

Code
par(op)

Notice how the comorbidity of obesity results in many more infected agents than the when the comorbidity isn’t present. Also, note how the plot_incidence() function output differs from that of the plot() function.

Note

If you want to drill further into this data, you can get the agents’ final states using the function get_agents_states().

Visualization Examples

Transmission Network using netplot and igraph

Code
suppressPackageStartupMessages(library(netplot))
suppressPackageStartupMessages(library(igraph))

model_sir <- ModelSIR(
  name              = "COVID-19",
  prevalence        = .01,
  transmission_rate = .5,
  recovery_rate     = .5
  ) |> 
    agents_smallworld(n = 500, k = 10, d = FALSE, p = .01) |>
    verbose_off() |>
    run(ndays = 50, seed = 1912)

## Transmission network
net <- get_transmissions(model_sir)

## Plot
x <- graph_from_edgelist(as.matrix(net[,2:3]) + 1)

nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex=TRUE)

Multiple Simulations

Code
model_sir <- ModelSIRCONN(
  name              = "COVID-19",
  prevalence        = 0.01,
  n                 = 1000,
  contact_rate      = 2,
  transmission_rate = 0.9,
  recovery_rate     = 0.1
  )

## Generating a saver
saver <- make_saver("total_hist", "reproductive")

## Running and printing
run_multiple(model_sir, ndays = 100, nsims = 50, saver = saver, nthread = 2)
Starting multiple runs (50)
_________________________________________________________________________
_________________________________________________________________________
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
 done.
Code
## Retrieving the results
mult_sim_results <- run_multiple_get_results(model_sir)

head(mult_sim_results$total_hist)
  sim_num date nviruses       state counts
1       1    0        1 Susceptible    990
2       1    0        1    Infected     10
3       1    0        1   Recovered      0
4       1    1        1 Susceptible    974
5       1    1        1    Infected     26
6       1    1        1   Recovered      0
Code
plot(mult_sim_results$reproductive)

Likelihood-Free Markhov Chain Monte Carlo (LFMCMC)

Often in network simulations, we don’t know the model parameters that will produce an accurate model. Likelihood-Free Markhov chain Monte Carlo (LFMCMC) runs a base model over a specified number of simulations, each time modifying the parameters used by the model to bring the model results closer to the observed data we’re trying to model. In epiworldR, the LFMCMC() function creates an LFMCMC object that can perform this calibration. In this example, we use LFMCMC to calibrate a COVID-19 SIR Model. For an example of using epiworldR and the LFMCMC() function specifically to calibrate a model on real-world data, see the epiworld-forecasts project.

Create a True COVID-19 Model

Assume that the true parameters of COVID-19 for a given population of 2,000 agents are as follows:

  • Initial Disease Prevalence: 0.01
  • Transmission Rate: 0.1
  • Recovery Rate: 1/7

We would represent that disease in epiworldR using the ModelSIR() and agents_smallworld() functions.

Code
library(epiworldR)
model_seed <- 122

true_covid_model <- ModelSIR(
    name = "COVID-19",
    prevalence = .01,
    transmission_rate = .1,
    recovery_rate = 1/7
    )

agents_smallworld(
  true_covid_model,
  n = 2000,
  k = 5,
  d = FALSE,
  p = 0.01
)

Running the true model for 50 days results in the following final distribution of agents across the three SIR states:

Code
verbose_off(true_covid_model)
run(true_covid_model, ndays = 50, seed = model_seed)
observed_data <- get_today_total(true_covid_model)
observed_data
Susceptible    Infected   Recovered 
       1865           0         135

For the rest of the example, we’ll assume we don’t know the true disease parameters, but that we have the observed_data (e.g., from public health records). We’ll use LFMCMC to recover the transmission and recovery rates from the true model and use the observed_data to check how close each simulation is to the true values.

Setup LFMCMC

Frist, set up a new SIR model for LFMCMC to use. Since we don’t know the true parameters, we’ll guess. It won’t matter what we choose for the recovery and transmission rates, as we’ll define the initial parameters for LFMCMC before running it.

Code
model_sir <- ModelSIR(
    name = "COVID-19",
    prevalence = .01,
    transmission_rate = .9, #TODO:?
    recovery_rate = .3
    )

agents_smallworld(
  model_sir,
  n = 2000,
  k = 5,
  d = FALSE,
  p = 0.01
)

verbose_off(model_sir)

Next, define the LFMCMC functions (described in more detail here). Since we are trying to recover the Transmission and Recovery rates, our simulation function will test a new set of those two parameters during each iteration of LFMCMC.

Code
# Define the simulation function
simulation_fun <- function(params, lfmcmc_obj) {
  set_param(model_sir, "Recovery rate", params[1])
  set_param(model_sir, "Transmission rate", params[2])
  run(model_sir, ndays = 50)
  res <- get_today_total(model_sir)
  return(res)
}

# Define the summary function
summary_fun <- function(data, lfmcmc_obj) {
  return(data)
}

# Define the proposal function
proposal_fun <- function(old_params, lfmcmc_obj) {
  res <- plogis(qlogis(old_params) + rnorm(length(old_params), sd = .1))
  return(res)
}

# Define the kernel function
kernel_fun <- function(
    simulated_stats, observed_stats, epsilon, lfmcmc_obj
    ) {

  diff <- ((simulated_stats - observed_stats)^2)^epsilon
  dnorm(sqrt(sum(diff)))
}

Finally, use the LFMCMC() function to create the LFMCMC object, add the LFMCMC functions defined above, and pass in the observed COVID-19 data.

Code
lfmcmc_model <- LFMCMC(model_sir) |>
  set_simulation_fun(simulation_fun) |>
  set_summary_fun(summary_fun) |>
  set_proposal_fun(proposal_fun) |>
  set_kernel_fun(kernel_fun) |>
  set_observed_data(observed_data)

Run LFMCMC

Before running LFMCMC, we need to pick the initial parameters (recovery rate = 0.3, transmission rate = 0.3) and choose an epsilon for the kernel function (epsilon = 1.0). We’ll run LFMCMC for 2000 iterations or samples (n_samples = 2000).

Code
# Set initial parameters
init_params <- c(0.3, 0.3)
epsilon <- 1.0
n_samples <- 2000

# Run the LFMCMC simulation
verbose_off(lfmcmc_model)
run_lfmcmc(
  lfmcmc = lfmcmc_model,
  params_init = init_params,
  n_samples = n_samples,
  epsilon = epsilon,
  seed = model_seed
)

Check Results

Print the LFMCMC object with a burn-in period of 1,500.

Code
set_stats_names(lfmcmc_model, get_states(model_sir))
set_params_names(lfmcmc_model, c("Recovery rate", "Transmission rate"))

print(lfmcmc_model, burnin = 1500)
___________________________________________

LIKELIHOOD-FREE MARKOV CHAIN MONTE CARLO

N Samples (total) : 2000
N Samples (after burn-in period) : 500
Elapsed t : 1.00s

Parameters:
  -Recovery rate     :  0.14 [ 0.13,  0.15] (initial :  0.30)
  -Transmission rate :  0.09 [ 0.09,  0.10] (initial :  0.30)

Statistics:
  -Susceptible :  1865.49 [ 1864.00,  1866.00] (Observed:  1865.00)
  -Infected    :     0.00 [    0.00,     0.00] (Observed:     0.00)
  -Recovered   :   134.51 [  134.00,   136.00] (Observed:   135.00)
___________________________________________

Out LFMCMC calibration was successful! The average LFMCMC Recovery rate was 0.14 (true value was 1/7 or 0.1428) and the average Transmission rate was 0.09 (true value was 0.1). The average number of Susceptible agents at Day 50 was 1,865.49 (observed value was 1865) and the average number of Recovered agents was 134.51 (observed value was 135).