Skip to contents

Toy Example

Source: vignettes/cowflu.Rmd
cowflu.Rmd

This vignette walks through the basic use of the model for a toy example.

Each model system is comprised of a number of regions, and each region contains a unique number of herds.

Consider 3 regions, each containing 3, 7, and 11, herds respectively.

##Setting model parameters

First, we set the model parameters:

pars <- cowflu:::cowflu_inputs(
  alpha = 0.2, #The rate of herd-to-herd infection contacts within each region.
  beta = 0.9,  #The transmission rate
  gamma = 0.7, #The incubation rate
  sigma = 0.9,  #The recovery rate
  asc_rate = 1, # The ascertainment rate
  dispersion = 1, # Dispersion parameter, when fitting
  inputs = cowflu:::cowflu_fixed_inputs(
    n_herds_per_region = c(3, 7, 11),
    p_region_export = c(.5, .5, .5),
    p_cow_export = c(0.2, 0.2, 0.2),
    n_cows_per_herd = c(rep(200, 3), rep(1000, 7), rep(3000, 11)),
    movement_matrix = cbind(c(.6, .2, .2), c(.2, .6, .2), c(.2, .2, .6)),
    time_test = 10000,
    start_herd = 5,
    start_count = 5))

alpha, beta, gamma, and sigma are epidemiological parameters. n_herds_per_region is a vector listing the number of herds in each region. p_region_export is a vector listing the probability that, each day, a herd within the respective region will export cattle to another herd. This does not yet inform where that herd will go, only whether or not an export is going to happen. p_cow_export is a vector listing the probability that, given a herd is exporting cattle, what is the probability that each cow in that herd will be included in the export. i.e. our toy example states that if a herd is exporting cattle, 20% of that herd will be exported on average, in every region. n_cows_per_herd is a vector listing exactly how many cows are in each herd of the model. movement_matrix details the probability of which region an exported herd will be sent to. Element (i, j) is the probability that an export from region i will be sent to region j. Thus, each row will sum to 1. time_test is the time point (in days) at which cattle will start being tested for flu. At this point, when a herd is exported from one region to another, a number of cows (see x below) will be sampled from the exported cows. If any of these sampled cows are in the Infected compartment, the export is cancelled. start_herd is the herd in which infected cows are seeded into at model initialisation. start_count is the number of infected cows to be seeded into herd start_herd at model initialisation.

Other model parameters are: n_test is the number of cows that are tested during inter-region exports from time_test onwards. This defaults to 30.

Running the model

Start by setting the seed, and the number of particles we wish to run

set.seed(1)
n_particles <- 1

Initialise the dust system:

sys <- dust2::dust_system_create(cowflu:::cows(), pars, n_particles = n_particles, dt = 0.25, seed = 42)
dust2::dust_system_set_state_initial(sys)

Run the model for the desired number of time-steps, in this instance, from 0 to 75:

s <- dust2::dust_system_simulate(sys, 0:75)

One of the model variables exported on the end of this s object is not part of the SEIR dynamics. We want to chop it off for now.

end_of_core_states <- (pars$n_herds + pars$n_regions)*5
s <- s[1:end_of_core_states,]

Now re-arrange this output into an array:

s1 <- array(s, c(pars$n_herds + pars$n_regions, 5, n_particles, 76))

Dim 1 is the herd, and note that we also export cumulative region totals too on the end.

Dim 2 is the number of epi compartments: S, E, I, and R, and outbreak recording.

Dim 3 is the number of particles run.

Dim 4 is time.

Note, for large models, we may not want (or be able) to export every herd’s dynamic profile. Instead, we can ask dust2 to just return the region totals like so:

Create i, a list of the region total indices.

i <- seq.int(pars$n_herds + 1, length.out = pars$n_regions)
i <- c(outer(i, (pars$n_herds + pars$n_regions) * (0:3), "+"))

And run the model, and allocate output into an array

sys <- dust2::dust_system_create(cowflu:::cows(), pars, n_particles = n_particles, dt = 0.25, seed = 42)
dust2::dust_system_set_state_initial(sys)
s <- dust2::dust_system_simulate(sys, 0:75, i)
s2 <- array(s, c(pars$n_regions, 4, n_particles, 76))

And check that these outputs are the same:

testthat::expect_equal(s1[22:24, 1:4, , , drop = FALSE], s2)