Part 3b: 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 part of the workshop, we’ll walk through an example.

Example Scenario: Calibrate a COVID-19 SIR Model

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.

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:

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.

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.

# 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.

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).

# 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.

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).