Commuter Model

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Background

This model is a stochastic SEIIaR (susceptible, exposed, infectious, infectious asymptomatic, recovered) metapopulation model. Each location has a local infection system, while the locations are connected by people who commute each day. The model differentiates between day and night. During the day you can infect/be infected in the location where you work, while during the night you can infect/be infected in the location where you live. It is the same commuters who travel back and forth each day. At the start of a day, all commuters are sent to their work location, where they mix for 12 hours. The commuters are then sent to their respective home locations, where they mix for 12 hours. The model is loosely based upon a published model by Engebretsen (2019).

Data required

library(ggplot2)
library(data.table)

seiiar is a dataset that contains the following variables for the entire population you want to model:

# no one in Norway is infected, and everyone is susceptible
spread::norway_seiiar_noinfected_2017

# 10 people in Oslo are infected, and everyone is susceptible
spread::norway_seiiar_oslo_2017

# no one in Norway is infected, and childhood vaccination data is used to
# estimate the number of "recovered" (i.e. non-susceptible) people for measles
spread::norway_seiiar_measles_noinfected_2017

# 10 people in Oslo is infected, and childhood vaccination data is used to
# estimate the number of "recovered" (i.e. non-susceptible) people for measles
spread::norway_seiiar_measles_oslo_2017

# we can take a closer look at Oslo
spread::norway_seiiar_measles_oslo_2017[location_code=="municip0301"]

commuters is a dataset that contains the following variables:

# we provide the number of municipal commuters in Norway in 2017
spread::norway_commuters_2017

r0 is the basic reproductive number.

latent_period is the average number of days from when a person is exposed until they are infectious.

asymptomatic_prob is the probability that an infectious person is asymptomatic.

asymptomatic_relative_infectiousness is the relative infectiousness of an asymptomatic infectious person when compared to a symptomatic infectious person. This value is between 0 and 1.

days_simulation the number of days you want to simulate.

N this is the number of simulations that will be run and then the results will be averaged. This should generally be set to 1, unless you are performing model fitting.

Example

We simulate one measles outbreak using the datasets spread::norway_seiiar_measles_oslo_2017 and spread::norway_commuters_2017:

set.seed(4)
d <- spread::commuter(
  seiiar=spread::norway_seiiar_measles_oslo_2017,
  commuters=spread::norway_commuters_2017,
  r0=14,
  latent_period = 8,
  infectious_period = 5,
  asymptomatic_prob=0,
  asymptomatic_relative_infectiousness=0,
  days_simulation=7*9,
  N=1
)

Inspecting the data from Oslo:

d[location_code=="municip0301"]

We can merge in information about counties, and create county level data:

d <- merge(d,fhidata::norway_locations_current, by.x="location_code",by.y="municip_code")
county <- d[,.(
  S=sum(S),
  E=sum(E),
  I=sum(I),
  Ia=sum(Ia),
  R=sum(R),
  incidence=sum(incidence),
  pop=sum(pop)
),
keyby=.(county_code,county_name,week,day,is_6pm)]
county[,county_name:=factor(county_name,levels=unique(fhidata::norway_locations_current[,c("county_code","county_name")]$county_name))]
county

We can produce a daily epicurve for each county:

p <- ggplot(county, aes(x=day, y=incidence))
p <- p + geom_col()
p <- p + facet_wrap(~county_name)
p <- p + scale_x_continuous("Day")
p

And we can produce weekly maps:

w <- county[,.(
  incidence_weekly = sum(incidence),
  pop = mean(pop)
), keyby=.(county_code, week)]

w[,weekly_incidence_per_10000 := 10000*incidence_weekly/pop]
w[,facet:=glue::glue("Week {week}",week=week)]

pd <- merge(
  w,
  fhidata::norway_map_counties, 
  by.x="county_code",
  by.y="location_code",
  allow.cartesian = T)

p <- ggplot(data=pd, mapping=aes( x = long, y = lat, group = group))
p <- p + geom_polygon(aes(fill=weekly_incidence_per_10000))
p <- p + facet_wrap(~facet)
p <- p + theme_void()
p <- p + coord_quickmap()
p


Try the spread package in your browser

Any scripts or data that you put into this service are public.

spread documentation built on Aug. 20, 2019, 5:13 p.m.