This repo contains a bunch of test code from various stages of making this work and most files in sub-directories can be ignored. The only relevant files now are those in the root directory and MiP-given/equilibrium-init-create.R, MiP-given/model_parameters.R, and MiP-given/MZ_multi_rates.rds for generating the model parameters.

Install the latest versions of odin.dust and mcstate from GitHub

remotes::install_github("mrc-ide/odin.dust", upgrade = FALSE)
remotes::install_github("mrc-ide/mcstate", upgrade = FALSE)

Run and plot particle trajectories by running the run-toy.R script. This follows these steps:

Prepare data set using mcstate::particle_filter_data. rate must be NULL here and initial_time will most likely be 0 (the first data value must be > 0, so the times in the original data provided have been incremented accordingly)

data_raw <- read.csv("casedata_monthly.csv",
                       stringsAsFactors = FALSE)
data <- mcstate::particle_filter_data(data_raw, time = "t", rate = NULL, initial_time = 0)

Define an index function for filtering the run state. The first list item, run should contain the portion of state needed for the likelihood calculation; the second list item, state should contain the portion of state for which particle history will be saved (see below)

index <- function(info) {
  list(run = c(Ih = info$index$Ih),
       state = c(Ih = info$index$Ih,
                 Sh = info$index$Sh))
}

Define a comparison function

compare <- function(state, observed, pars = NULL) {
    Ih <- state[1, ] # as defined by the above index
    dbinom(x = observed$positive,
           size = observed$tested,
           prob = Ih,
           log = TRUE)
}

A schedule for running the stochastic updates

stochastic_schedule <- seq(from = 60, by = 30, to = 1830)

Run the model and get a likelihood

model <- odin.dust::odin_dust("toyodinmodel.R")
n_particles <- 100
p <- mcstate::particle_filter$new(data, model, n_particles, compare,
                       index = index, seed = 1L,
                       stochastic_schedule = stochastic_schedule)
pars <- list(init_Ih = 0.8,
             init_Sv = 100,
             init_Iv = 1,
             nrates = 15)
lik <- p$run(pars)

To plot particle trajectories, run the model with save_history = TRUE:

lik <- p$run(pars, save_history = TRUE)
history <- p$history()
matplot(data_raw$t, t(history[1, , -1]), type = "l",
        xlab = "Time", ylab = "State",
        col = "#ff000022", lty = 1, ylim = range(history))
matlines(data_raw$t, t(history[2, , -1]), col = "#0000ff22", lty = 1)

Run and plot particle trajectories by running the run.R script. This follows exactly the same steps as above. betaa_td in this model isn't actually being updated via a stochastic process, but just incremented every time step as per the example you gave us. This can be edited in the model code mipodinmodel.R (see the toy model syntax in toyodinmodel.R as a reference). The parameters for the model have been generated using the exact same scripts you provided us (now found at MiP-given/model_parameters.R and MiP-given/equilibrium-init-create.R)