wfmackey
commented
4 years ago Playing around with this now. @HughParsonage can you tell me what is happening with initial case numbers for a given by_state = date? Looks very different to the covid19.model.sa2/inst/extdata/time_series_cases.fst file.
library(tidyverse)
options(dplyr.summarise.inform = FALSE)
library(purrr)
library(data.table)
#>
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:dplyr':
#>
#> between, first, last
#> The following object is masked from 'package:purrr':
#>
#> transpose
library(patchwork)
library(covid19.model.sa2)
library(fst)
library(glue)
#>
#> Attaching package: 'glue'
#> The following object is masked from 'package:dplyr':
#>
#> collapse
library(lubridate)
#>
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:data.table':
#>
#> hour, isoweek, mday, minute, month, quarter, second, wday, week,
#> yday, year
#> The following objects are masked from 'package:base':
#>
#> date, intersect, setdiff, union
library(janitor)
#>
#> Attaching package: 'janitor'
#> The following objects are masked from 'package:stats':
#>
#> chisq.test, fisher.test
run_days <- 30
repeat_times <- 1
repeat_simulation <- function(.name,
.repeat = repeat_times,
.days = run_days,
PolicyPars = set_policypars(),
EpiPars = set_epipars(),
...) {
return_simulation <- function(x,
to_return = .to_return,
.PolicyPars = PolicyPars,
.EpiPars = EpiPars) {
rep_each <- function(x, len) {
rep(x, each = len / length(x), length.out = len)
}
d <-
simulate_sa2(days_to_simulate = .days,
returner = 3,
PolicyPars = .PolicyPars,
EpiPars = .EpiPars)
# Status
Status12 <- d[["Status12"]]
DT <- data.table(N = Status12)
DT[, day := rep_each(1:.days, .N)]
DT[, state := rep_each(1:9, .N), by = .(day)]
DT[, status := rep_each(-2:3, .N), by = .(day, state)]
DT[, isolated := rep_len(0:1, .N)]
DT[, runid := x]
# New infections
NewInfectionsByState <- d[["NewInfectionsByState"]]
DT2 <- data.table(new_infections = NewInfectionsByState)
DT2[, day := rep_each(1:.days, .N)]
DT2[, state := rep_each(1:9, .N), by = .(day)]
return(left_join(DT, DT2))
}
d <- map_dfr(1:.repeat, return_simulation)
if (!dir.exists("output")) dir.create("output")
write_fst(d, glue("output/{.name}.fst"))
}
plot_status <- function(data, keep_all_states = FALSE) {
if (!keep_all_states) {
data <- data %>%
filter(state <= 6 & state != 5)
}
data %>%
filter(status >= 1) %>%
group_by(runid, state, day) %>%
summarise(N = sum(N)) %>%
ggplot(aes(day, N)) +
geom_line() +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous(limits = c(0, NA)) +
facet_wrap(vars(state), scales = "free_y", ncol = 1) +
theme(legend.position = "top")
}
plot_new_cases <- function(data, keep_all_states = FALSE) {
if (!keep_all_states) {
data <- data %>%
filter(state <= 6 & state != 5)
}
data %>%
group_by(runid, state, day) %>%
filter(row_number() == 1) %>%
ggplot(aes(day, new_infections)) +
geom_line() +
facet_wrap(vars(state), scales = "free_y", ncol = 1)
}
# Looking at early days
from_date <- ymd("2020-03-05")
keep_states <- c("NSW", "Vic")
early_days <- repeat_simulation("early_days",
.repeat = 1,
.days = 21,
by_state = from_date,
EpiPars = set_epipars(a_workplace_rate = .5,
q_workplace = .01,
q_household = 0.02,
a_household_rate = 0.15, # def 0.15
q_supermarket = 0.002),
PolicyPars = set_policypars(
supermarkets_open = TRUE,
max_persons_per_supermarket = 100,
schools_open = TRUE,
only_Year12 = FALSE,
school_days_per_wk = 5,
do_contact_tracing = FALSE,
workplaces_open = TRUE,
workplace_size_max = 100
))
#> Joining, by = c("day", "state")
# Tidy early_days
early_days_data <- early_days %>%
mutate(state = strayr::strayr(state),
infected = status >= 1) %>%
filter(infected) %>%
group_by(runid, state, day) %>%
summarise(cases = sum(N),
new_cases = max(new_infections)) %>%
mutate(date = from_date + (day - 1)) %>%
filter(state %in% keep_states)
# Get actual case numbers
actual <- read_fst("../covid19.model.sa2/inst/extdata/time_series_cases.fst") %>%
clean_names() %>%
select(-total) %>%
pivot_longer(-date, names_to = "state", values_to = "cases") %>%
mutate(state = strayr::strayr(state)) %>%
filter(state %in% keep_states) %>%
group_by(state) %>%
mutate(new_cases = cases - lag(cases))
# Plot
early_days_data %>%
ggplot(aes(x = date, y = cases)) +
geom_line(aes(group = runid),
alpha = .8, colour = "red") +
geom_line(data = actual, colour = "black") +
facet_wrap(vars(state)) +
scale_x_date(limits = c(from_date, from_date + 21)) +
scale_y_continuous() +
theme_minimal()
#> Warning: Removed 97 row(s) containing missing values (geom_path).
Created on 2020-05-26 by the reprex package (v0.3.0)
HughParsonage
This will be shown at the end of Chapter 3.