2.4 Desiderata for epidemiological models

library(here)
source(here("setup.R"))

here() starts at /Users/stefan/workspace/work/phd/thesis

w_df <- read_csv(here("data/processed/generation_time.csv"))
w <- w_df$w

estimate_R <- function(I, w) {
    I_ext <- c(rep(0, length(w) - 1), I)
    I / stats::filter(I_ext, c(0, w), sides = 1, method = "convolution")[-(1:(length(w) - 1))]
}

estimate_rho <- function(I) {
    I / lag(I, 7)
}

fct_events <- ordered(c("outbreak", "christmas"))

stopifnot(length(estimate_R(rpois(100, 10), w)) == 100)

rki <- read_csv(here("data/processed/rki_county.csv"))


rki %>%
    group_by(date) %>%
    summarize(cases = sum(cases)) %>%
    mutate(cases_7 = rollmean(cases, k = 7, fill = NA)) %>%
    mutate(R = estimate_R(cases, w)) %>%
    mutate(R_7 = estimate_R(cases_7, w)) %>%
    mutate(rho = estimate_rho(cases)) %>%
    mutate(rho_7 = estimate_rho(cases_7)) %>%
    filter(date >= "2020-04-01") %>%
    select(date, R:rho_7) %>%
    pivot_longer(-date) %>%
    mutate(type = ifelse(str_starts(name, "R"), "$\\hat R$", "$\\hat \\rho^7$")) %>%
    mutate(type = factor(type, levels=c( "$\\hat \\rho^7$", "$\\hat R$"))) %>%
    mutate(timescale = ifelse(str_ends(name, "_7"), "7-day avg.", "daily")) %>%
    ggplot(aes(date, value, color = timescale, alpha = timescale, group = timescale)) +
    geom_line() +
    facet_wrap(~type, scales = "free_y", nrow = 2) +
    scale_color_manual(values = c("daily" = "black", "7-day avg." = pal_npg()(1))) +
    scale_alpha_manual(values = c("daily" = .4, "7-day avg." = 1)) +
    coord_cartesian(ylim = c(0, 3.2)) +
    scale_x_date(date_minor_breaks = "1 month", limits = c(ymd("2020-04-01"), ymd("2021-03-01")), date_breaks = "3 months", date_labels = "%b %Y") +
    theme(legend.position = "bottom") +
    geom_mark_rect(
        aes(
            x = date, y = value,
            filter = (date <= ymd("2020-07-01") & date >= ymd("2020-06-01")), group = type, linetype = fct_events[1]
        ),
        inherit.aes = FALSE
    ) +
    geom_mark_rect(
        aes(
            x = date, y = value,
            filter = (date >= ymd("2020-12-07") & date <= ymd("2021-01-14")),
            group = type, linetype = fct_events[2]
        ),
        inherit.aes = FALSE
    ) +
    scale_linetype_manual(values = c("outbreak" = "dashed", "christmas" = "dotted")) +
    labs(x = "", y = "", color = "incidences", alpha = "incidences", linetype = "")
# geom_magnify( # Toennies
#    from = list(ymd("2020-06-01"), ymd("2020-07-01"), .6, 2.2),
#    to = list(ymd("2020-08-01"), ymd("2020-10-14"), 2.1, 3),
#    axes = "x"
# )
# geom_magnify( # christmas
#    from = list(ymd("2020-12-07"), ymd("2021-01-14"), .5, 1.5),
#    to = list(ymd("2021-03-01"), ymd("2021-07-01"), 2.2, 3.9),
#    axes = "x"
# )

ggsave_tikz(here("tikz/rho_and_R_naive.tex"))

Warning message:

“Removed 1590 rows containing missing values or values outside the scale range

(`geom_line()`).”

Warning message:

“Removed 1590 rows containing missing values or values outside the scale range

(`geom_line()`).”

pdf: 2

rho <- seq(0.01, 1.1, length.out = 1001)

rho_to_R <- function(rho, w) {
    1 / sum(rho^(-seq_along(w)) * w)
}
R_c <- sapply(rho, function(rho) rho_to_R(rho, w_df$w))
EW <- sum(seq_along(w) * w)

tibble(rho = rho, R_c = R_c, rho_7 = rho^7, rho_EW = rho^EW) %>%
    filter(R_c <= 1.2, R_c >= 0.8) %>%
    summarize(
        max_diff_rho_7 = max(abs(R_c - rho_7)),
        max_diff_rho_EW = max(abs(R_c - rho_EW))
    )

A tibble: 1 × 2

max_diff_rho_7	max_diff_rho_EW
<dbl>	<dbl>
0.05675741	0.002564586

5.61111111111111