-
Notifications
You must be signed in to change notification settings - Fork 10
Household SEIR
In this example, we will demonstrate how to define a contact pattern. The package currently provide a random-mixing contact pattern, and a random contact network generated by a configuration model. We can define new contact patterns by subclassing an R6 class named Contact.
Here we define a household contact. Each household holds a random number of individuals, which is generated using an R function. This function is passed to the constructor.
There are two crucial methods that need to be defined: $build()
and $contact(time, agent)
.
-
$build()
is called when the contact pattern is attached to a population, and it is the time when the contact of each individual is defined. When this method is called, theprivate$population
field is defined and holds an external pointer to the population. Here we randomly create households, and put the agents in the population into these households.- Households are a list of external pointers pointing to the ants in the household, and is stored in
$private$households
- For each agent, all agents in the same households are the contacts.
- We associate the households with the agents' id, so that we can quickly lookup the households to which each agent belongs. This association is stored in
private$agents
.
- Households are a list of external pointers pointing to the ants in the household, and is stored in
-
$contact(time, agent)
is called when an agent requests a contact at the given time. It should return a list of agents to contact. In this example, the contacts are the agents in the same household.
Household = R6::R6Class(
"Household",
inherit = Contact,
public = list(
initialize = function(rng) {
super$initialize()
private$rng = rng
},
build = function() {
n = getSize(private$population)
agents = rep(NA, n)
k = 1
h = 1
ok = TRUE
while (k <= n) {
m = private$rng()
hh = list()
for (j in 1:m) {
hh[[j]] = getAgent(private$population, k)
private$agents[k] = h
k = k + 1
if (k > n) break
}
private$households[[h]] = hh
h = h + 1
}
},
contact = function(time, agent) {
id = getID(agent)
private$households[[private$agents[[id]]]]
}
),
private = list(
rng = NULL,
households = NULL,
agents = NULL
)
)
As population can have multiple contact patterns attached. In this example, agents contact other agents in the same household at an per-contact rate beta_hh, and also randomly contact an agent at a per-capita rate beta_rm.
In addition, we use gamma distributed latent and infectious periods in this example. The code is almost identical to the SEIR model
gamma = newGammaWaitingTime(5, 1)
beta.rm = 0.05
beta.hh = 0.1
sigma= newGammaWaitingTime(5, 1)
n = 10000
I0 = 10
sim = Simulation$new(as.list(c(rep("I", I0), rep("S", n-I0))))
sim$addLogger(newCounter("S", "S"))
sim$addLogger(newCounter("E", "E"))
sim$addLogger(newCounter("I", "I"))
sim$addLogger(newCounter("R", "R"))
m = newRandomMixing()
sim$addContact(m)
h = Household$new(function() { rpois(1, 7) + 1})
sim$addContact(h)
sim$addTransition("E"->"I", sigma)
sim$addTransition("I"->"R", gamma)
sim$addTransition("I" + "S" -> "I" + "E" ~ m, beta.rm)
sim$addTransition("I" + "S" -> "I" + "E" ~ h, beta.hh)
sim$run(0:300)