Simulates cases by date of report from input Rt and initial case seed

simulate-cases

library(EpiNow)
library(EpiSoon)
library(data.table)

#1. Input a vector of Rt estimates, and vector to initialise cases, and a generation time
initial_cases <- 100
initial_date <- as.Date("2020-03-01")
rts <- c(rep(2, 20), (2 - 1:15 * 0.1), rep(0.5, 10))

## Simulating cases
cases <- data.table::data.table(
  date = initial_date,
  cases = initial_cases)

## Structuring rs
rts <- data.table::data.table(
  date = seq(initial_date, initial_date + lubridate::days(length(rts) - 1), by = "days"),
  rt = rts
)

generation_interval <- rowMeans(EpiNow::covid_generation_times)
median_interval <- sum(!(cumsum(generation_interval) > 0.5)) + 1
#2 Return a dataframe of cases by date of infection
simulated_cases <- EpiSoon::predict_cases(cases = cases,
                                          rt = rts,
                                          serial_interval = generation_interval,
                                          rdist = rpois) 

#3 Input cases date of infection, report delay (with uncertainty) and incubation period (with uncertainty)
## Define a single report delay distribution
delay_def <- EpiNow::lognorm_dist_def(mean = 5,
                                      mean_sd = 1,
                                      sd = 3,
                                      sd_sd = 1,
                                      max_value = 30,
                                      samples = 1)

## Define a single incubation period
incubation_def <- EpiNow::lognorm_dist_def(mean = EpiNow::covid_incubation_period[1, ]$mean,
                                           mean_sd = EpiNow::covid_incubation_period[1, ]$mean_sd,
                                           sd = EpiNow::covid_incubation_period[1, ]$sd,
                                           sd_sd = EpiNow::covid_incubation_period[1, ]$sd_sd,
                                           max_value = 30, samples = 1)

## Mapping with a weekly reporting effect
report_weekly <- adjust_infection_to_report(
                      simulated_cases, delay_def, incubation_def,
                      reporting_effect = c(1.1, rep(1, 4), 0.95, 0.95))


## Truncate initial data be length of the median generation interval
report_weekly[date >= min(date) + lubridate::days(median_interval)]
s#4. Return cases by date of report (optionally also return cases by date of onset?)

print(report_weekly)

Now integrated here: https://epiforecasts.io/EpiNow/reference/simulate_cases.html