Last active
August 21, 2020 17:26
-
-
Save seabbs/07035fa4019a0c0117e28f9037593f28 to your computer and use it in GitHub Desktop.
Evaluates Rt over two fixed time points with Rt assumed to vary based on a random walk in-between the two time points and after a period of time.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library(EpiNow2) | |
| library(data.table) | |
| # Get example case counts | |
| reported_cases <- EpiNow2::example_confirmed[1:60] | |
| # Add a dummy breakpoint (used only when optionally estimating breakpoints) | |
| reported_cases <- reported_cases[, breakpoint := data.table::fifelse((date >= as.Date("2020-03-16") & date <= as.Date("2020-03-29")) | date >= as.Date("2020-04-11"), | |
| 1, 0)] | |
| # Set up example generation time | |
| generation_time <- list(mean = EpiNow2::covid_generation_times[1, ]$mean, | |
| mean_sd = EpiNow2::covid_generation_times[1, ]$mean_sd, | |
| sd = EpiNow2::covid_generation_times[1, ]$sd, | |
| sd_sd = EpiNow2::covid_generation_times[1, ]$sd_sd, | |
| max = 30) | |
| # Set delays between infection and case report | |
| # (any number of delays can be specifed here) | |
| incubation_period <- list(mean = EpiNow2::covid_incubation_period[1, ]$mean, | |
| mean_sd = EpiNow2::covid_incubation_period[1, ]$mean_sd, | |
| sd = EpiNow2::covid_incubation_period[1, ]$sd, | |
| sd_sd = EpiNow2::covid_incubation_period[1, ]$sd_sd, | |
| max = 30) | |
| reporting_delay <- list(mean = log(5), | |
| mean_sd = log(2), | |
| sd = log(2), | |
| sd_sd = log(1.5), | |
| max = 30) | |
| # Run model with default settings | |
| def <- estimate_infections(reported_cases, family = "negbin", | |
| generation_time = generation_time, | |
| delays = list(incubation_period, reporting_delay), | |
| samples = 1000, warmup = 200, cores = 4, | |
| horizon = 0, chains = 4, estimate_rt = TRUE, | |
| verbose = TRUE, return_fit = TRUE) | |
| # Plot output | |
| plots <- report_plots(summarised_estimates = def$summarised, | |
| reported = reported_cases) | |
| plots$summary | |
| # Run model with breakpoints but otherwise static Rt | |
| # This formulation may increase the apparent effect of the breakpoint but needs to be tested using | |
| # model fit criteria (i.e LFO). | |
| fbkp <- estimate_infections(reported_cases, family = "negbin", | |
| generation_time = generation_time, | |
| delays = list(incubation_period, reporting_delay), | |
| samples = 1000, warmup = 200, cores = 4, | |
| chains = 4, estimate_breakpoints = TRUE, fixed = TRUE, horizon = 0, | |
| verbose = TRUE, return_fit = TRUE) | |
| # Plot output | |
| plots <- report_plots(summarised_estimates = fbkp$summarised, | |
| reported = reported_cases) | |
| plots$summary | |
| # Pull out breakpoint summary | |
| fbkp$summarised[variable == "breakpoints"] | |
| # Pull out R summary | |
| fbkp$summarised[variable == "R"] |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Breakpoints
Rt estimates: