This vignette describes how changes in the reproduction number can be modelled with just case data and with a discrete serial interval.
Data from the 1918 Flu pandemic is used as an example. This data is provided in the R package @EpiEstim.
library(epidemia) library(EpiEstim) library(plotly) data("Flu1918") options(mc.cores = parallel::detectCores()) print(Flu1918)
## $incidence
## [1] 5 1 6 15 2 3 8 7 2 15 4 17 4 10 31 11 13 36 13
## [20] 33 17 15 32 27 70 58 32 69 54 80 405 192 243 204 280 229 304 265
## [39] 196 372 158 222 141 172 553 148 95 144 85 143 87 73 70 62 116 44 38
## [58] 60 45 60 27 51 34 22 16 11 18 11 10 8 13 3 3 6 6 13
## [77] 5 6 6 5 5 1 2 2 3 8 4 1 2 3 1 0
##
## $si_distr
## [1] 0.000 0.233 0.359 0.198 0.103 0.053 0.027 0.014 0.007 0.003 0.002 0.001flu <- Flu1918
First, we form the data argument to epim. Recall that this must store all observations and covariates needed to fit the model. In our case, this is easy because there are no covariates, and there is just one observation type - case data.
data <- data.frame( country = "A", cases = c(NA, flu$incidence), # pad start fludate = as.Date("1918-01-01") + seq(0, along.with=c(NA,flu$incidence)) ) # needed for weekly random walk data$week <- format(data$fludate, "%V")
A week column as been added to data because we will parameterise the reproduction number as a weekly random walk. This is done as follows.
# model the rep number as weekly random walk (no covariates) rt <- epirt( formula = R(country, fludate) ~ 0 + rw(time=week, prior_scale=0.1), r0=3 # prior expected reproduction number )
We will also need to specify a model for how the observations (daily cases) are observed from the latent infection series. This is done through the epiobs function.
# model observed cases as a proportion of infections constant over time cases <- epiobs( formula = cases(country, fludate) ~ 1, prior_intercept = rstanarm::normal(location=1, scale=0.01), link = "identity", i2o = rep(.25,4) )
In this case, we assume that everybody infected is recorded as a case, and that this happens with equal probability in the four days subsequent to infection.
We will use the NUTs sampler to fit the model. It is useful to explicitly control the parameters for sampling.
# use NUTs with given parameters sampling_args <- list( iter=1000, control=list(adapt_delta=0.95,max_treedepth=15), seed=12345 )
We are ready to fit the model. This is done through a call to epim.
# fit the model with various model params fm <- epim( rt = rt, obs = list(cases), data = data, pops = data.frame(country="A", pop=1e6), si = flu$si_distr, seed_days = 6, prior_tau = rstanarm::exponential(rate=4), algorithm = "sampling", sampling_args = sampling_args )
The resulting fit can be examined using various plotting functions provided by epidemia.
subplot( plot_rt(fm, plotly=T, smooth=7), plot_infections(fm, plotly=T), plot_obs(fm, type = "cases", plotly=T), plot_infectious(fm, plotly=T), shareX = T, nrows=2, margin=0.04 )