project_seir: Make a projection with a SEIR fit

In seananderson/covidseir: Bayesian SEIR Modelling for Multivariate COVID-19 Case Data

Description

Predict from a fit_seir() object, possibly with a future prediction. By default, the projection uses the estimated f values (fraction of normal contacts encountered for those physical distancing). The function also includes functionality to specify a vector of fixed f values starting on a given future date.

Usage

project_seir(
  obj,
  forecast_days = 0,
  f_fixed_start = NULL,
  f_fixed = NULL,
  f_multi = NULL,
  f_multi_seg = NULL,
  iter = seq_along(obj$post$R0),
  return_states = FALSE,
  imported_cases = 0,
  imported_window = 1,
  parallel = TRUE,
  X = obj$stan_data$X,
  ...
)

Arguments

obj	Output from fit_seir().
forecast_days	Number of projection days.
f_fixed_start	Optional day to start changing f. Must be set if f_fixed or f_multi is set. This is the day number to start, not the number of days into the projection. E.g., if 100 days were fitted to and f_fixed should start 5 days into the future, then f_fixed_start = 105.
f_fixed	An optional vector of fixed f values for the projection. Should be length forecast_days - (f_fixed_start - nrow(daily_cases) - 1). I.e. one value per day after f_fixed_start day.
f_multi	Multiplicative vector of f values. Same structure as f_fixed.
f_multi_seg	Which f to use for f_multi.
iter	MCMC iterations to include. Defaults to all.
return_states	Logical for whether to return the ODE states.
imported_cases	Number of cases to import starting on first projection day.
imported_window	Number of days over which to distribute imported cases.
parallel	Use parallel processing via future and furrr?
X	An optional model matrix that acts additively on log expected cases.
...	Other arguments to pass to rstan::sampling().

Details

Set a future::plan() and this function will operate in parallel across MCMC iterations using future and furrr.

Value

A data frame:

day

Day

data_type

Data-type column from the case data

mu

Expected number of cases

y_rep

Posterior predictive replicate observation

phi

Posterior draw of phi, the NB2 dispersion parameter, if included

.iteration

The MCMC iteration

Examples

cases <- c(
  0, 0, 1, 3, 1, 8, 0, 6, 5, 0, 7, 7, 18, 9, 22, 38, 53, 45, 40,
  77, 76, 48, 67, 78, 42, 66, 67, 92, 16, 70, 43, 53, 55, 53, 29,
  26, 37, 25, 45, 34, 40, 35
)
# Example fixed sample fractions:
s1 <- c(rep(0.1, 13), rep(0.2, length(cases) - 13))

# To use parallel processing with multiple chains:
# options(mc.cores = parallel::detectCores() / 2)

# Using the MAP estimate for speed in this example:
m <- fit_seir(
  cases,
  iter = 100,
  i0_prior = c(log(8), 0.2),
  samp_frac_fixed = s1,
  fit_type = "optimizing"
)
print(m)

p <- project_seir(m)
p

obs_dat <- data.frame(day = seq_along(cases), value = cases)

library(magrittr) # for %>%
tidy_seir(p) %>% plot_projection(obs_dat = obs_dat)

tidy_seir(p) %>%
  plot_residuals(obs_dat = obs_dat, obj = m)

# for parallel processing (optional)
# future::plan(future::multisession)
p <- project_seir(m,
  forecast_days = 100,
  f_fixed_start = 53,
  f_fixed = c(rep(0.7, 60), rep(0.2, 30)),
  iter = 1:25
)
p
tidy_seir(p) %>% plot_projection(obs_dat = obs_dat)

# fake example to show optional Rt colouring:
proj <- tidy_seir(p)
plot_projection(proj, obs_dat = obs_dat, Rt = rlnorm(nrow(proj)), col = "#00000050")


# Get threshold for increase:
# future::plan(future::multisession) # for parallel processing (optional)
# (only using 40 iterations for a fast example)
thresh <- get_threshold(m, iter = 1:40, show_plot = TRUE)
mean(thresh)

# Get doubling time (prevalence is declining in this example)
get_doubling_time(m, iter = 1:40, show_plot = TRUE)

states_with_Rt <- get_rt(m)
states_with_Rt

library(ggplot2)
ggplot(states_with_Rt, aes(time, Rt, group = .iteration)) +
  geom_line(alpha = 0.5)

seananderson/covidseir documentation built on June 3, 2021, 6:49 a.m.