seabbs/change_with_rw_padding

## change_with_rw_padding
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"]
	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"]