index.Rmd

---
title: "Orange County, CA COVID Situation Report"
site: workflowr::wflow_site
output:
  workflowr::wflow_html:
    toc: false
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F, fig.align = "center", fig.width=16,  message=F)
#knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
knitr::opts_chunk$set(autodep = TRUE)
library(lubridate)
library(here)
library(fs)
library(tidyverse)
library(tidybayes)
library(scales)
library(patchwork)
library(glue)
library(zoo)
theme_set(theme_bw(base_size =22)) 
```

```{r reading-model-objects, cache=FALSE}
full_county_popsize <- 3175692

model_name <- "SEIeIpRD"
loc_name <- "oc"

source(here("code", model_name, "helper_functions.R"))
source(here("code", model_name, "plot_functions.R"))

oc_data <- read.csv(here("data", "oc_data.csv"))
#if you want to run the analysis, change to TRUE. Otherwise the site will display the most recent results
run_analysis <- TRUE

if (run_analysis==TRUE){
library(rstan)
  
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

# Read OCHCA Data ---------------------------------------------------------
# oc_data <- read.csv(here("data", "oc_data.csv"))

oc_data$posted_date <- as.Date(oc_data$posted_date)

old_folder_name = "2020-10-18_2020-11-22"

#read in old results, this is used to calculate parameters of initial condition distributions
oc_post_old <- read_rds(here("code", model_name, loc_name, old_folder_name, "oc_post.rds"))
oc_model_objects_old <- read_rds(here("code", model_name, loc_name, old_folder_name, "model_objects.rds"))

#set start and end date of analysis
start_date <- as.Date("2020-10-25")
end_date <- as.Date("2020-11-29")

# Setup Model Objects -----------------------------------------------------

#acquire posterior samples for fraction infected, fraction infectious, fraction early infectious 
#nearest day to when current analysis starts
oc_post_old_wide <- oc_post_old %>%
  spread_draws(epi_curves[i]) %>%
  ungroup() %>%
  select(i, epi_curves, .draw) %>%
  mutate(g = (i - 1) %% oc_model_objects_old$n_curves + 1) %>%
  mutate(g_text = oc_model_objects_old$epi_curve_names[g]) %>%
  mutate(g_text = fct_reorder(g_text, g)) %>%
  mutate(t = oc_model_objects_old$dates_pp[ceiling(i / oc_model_objects_old$n_curves)]) %>%
  select(-i, - g) %>%
  pivot_wider(names_from = g_text, values_from = epi_curves)

target_date <- unique(oc_post_old_wide$t)[which.min(abs(unique(oc_post_old_wide$t) - start_date))]

for_estimating <- oc_post_old_wide %>%
  filter(t == target_date) %>%
  transmute(S=S,
            frac_carrs = 1 - S,
            frac_carrs_infec = (IE + IP) / (E + IE + IP),
            frac_infec_early = IE /(IE + IP))

#purpose: provide method of moments estimators for parameters of a beta distribution
#param: x: draws from a beta distribution
#output: method of moment estimates for alpha and beta of a beta distribution
mom_beta <- function(x) {
  x_bar <- mean(x)
  v_bar <- var(x)

  c(alpha = x_bar * (x_bar * (1 - x_bar) / v_bar - 1),
    beta = (1 - x_bar) * (x_bar * (1 - x_bar) / v_bar - 1))
}

#adjusted method of moments estimator for sensitivity analyses (mean twice as large)
mom_beta3 <- function(x) {
  x_bar <- 2*mean(x)
  v_bar <- var(x)
  
  c(alpha = x_bar * (x_bar * (1 - x_bar) / v_bar - 1),
    beta = (1 - x_bar) * (x_bar * (1 - x_bar) / v_bar - 1))
}

#adjusted method of moments estiamtor for sensitivity analyses (mean 1/2 as large)
mom_beta4 <- function(x) {
  x_bar <- .5*mean(x)
  v_bar <- var(x)
  
  c(alpha = x_bar * (x_bar * (1 - x_bar) / v_bar - 1),
    beta = (1 - x_bar) * (x_bar * (1 - x_bar) / v_bar - 1))
}

frac_carrs_prior <- mom_beta(for_estimating$frac_carrs)
frac_carrs_infec_prior <- mom_beta(for_estimating$frac_carrs_infec)
frac_infec_early_prior <- mom_beta(for_estimating$frac_infec_early)

init_func <- function() list(frac_carrs = frac_carrs_prior[1] / sum(frac_carrs_prior),
                             frac_carrs_infec = frac_carrs_infec_prior[1] / sum(frac_carrs_infec_prior),
                             frac_infec_early = frac_infec_early_prior[1] / sum(frac_infec_early_prior))

control_list <- list(adapt_delta = 0.99,
                     max_treedepth = 20)
# control_list <- list()

# update population size with median of the posterior for fraction susceptibles day before analysis start
new_pop_size <- as.integer(median(for_estimating$S) * 3175692L)

test <-   oc_data%>%
  lump_ochca_covid(time_interval_in_days = 3,
                           first_day = start_date,
                           last_day = end_date)

#create model objects to input into stan, see helper_functions.R for details
model_objects <- create_model_objects(ochca_covid = oc_data,
                                      first_day = start_date,
                                      last_day = end_date,
                                      time_interval_in_days = 3,
                                      forecast_in_days = 48,
                                      future_tests_per_day = 700,
                                      popsize = new_pop_size,
                                      frac_carrs_beta = frac_carrs_prior,
                                      frac_carrs_infec_beta = frac_carrs_infec_prior,
                                      frac_infec_early_beta = frac_infec_early_prior)

#create a model objects for smapling from the prior distributions, used in figure creation
model_objects_priors_only <- model_objects
model_objects_priors_only$priors_only <- T
model_objects_priors_only$forecast_in_days <- 0

#create new model objects for the 4 sensitivity analyses shown in the report
model_objects_sensitivity1 <- model_objects
model_objects_sensitivity1_priors_only <- model_objects_priors_only

model_objects_sensitivity1$log_R0_normal <- c(log(2.5), .5)
model_objects_sensitivity1_priors_only$log_R0_normal <- c(log(2.5), .5)

model_objects_sensitivity2 <- model_objects
model_objects_sensitivity2_priors_only <- model_objects_priors_only

model_objects_sensitivity2$log_R0_normal <- c(log(.53), .78)
model_objects_sensitivity2_priors_only$log_R0_normal <- c(log(.53), .78)

model_objects_sensitivity3 <- model_objects
model_objects_sensitivity3_priors_only <- model_objects_priors_only

model_objects_sensitivity3$frac_carrs_beta <- mom_beta3(for_estimating$frac_carrs)
model_objects_sensitivity3_priors_only$frac_carrs_beta <- mom_beta3(for_estimating$frac_carrs)

model_objects_sensitivity4 <- model_objects
model_objects_sensitivity4_priors_only <- model_objects_priors_only

model_objects_sensitivity4$frac_carrs_beta <- mom_beta4(for_estimating$frac_carrs)
model_objects_sensitivity4_priors_only$frac_carrs_beta <- mom_beta4(for_estimating$frac_carrs)


# Setup Results Folder ----------------------------------------------------
folder_name <- str_c(model_objects$actual_first_day, model_objects$actual_last_day, sep = "_")
folder_loc <- here("code", model_name, "oc", folder_name)
dir_create(folder_loc)

write_rds(model_objects, file.path(folder_loc, "model_objects.rds"))
write_rds(model_objects_sensitivity1, file.path(folder_loc, "model_objects_sensitivity1.rds"))
write_rds(model_objects_sensitivity2, file.path(folder_loc, "model_objects_sensitivity2.rds"))
write_rds(model_objects_sensitivity3, file.path(folder_loc, "model_objects_sensitivity3.rds"))
write_rds(model_objects_sensitivity4, file.path(folder_loc, "model_objects_sensitivity4.rds"))

# Sample Posterior --------------------------------------------------------
oc_post <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                data = model_objects,
                init = init_func,
                seed = 0,
                chains = 4,
                control = control_list)

write_rds(oc_post, file.path(folder_loc, "oc_post.rds"))


# Sample Prior ------------------------------------------------------------
oc_prior <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                 data = model_objects_priors_only,
                 init = init_func,
                 seed = 0,
                 chains = 4,
                 control = control_list)

write_rds(oc_prior, file.path(folder_loc, "oc_prior.rds"))


# Sample Posterior Sensitivity 1 ------------------------------------------
oc_post_sensitivity1 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                             data = model_objects_sensitivity1,
                             init = init_func,
                             seed = 0,
                             chains = 4,
                             control = control_list)

write_rds(oc_post_sensitivity1, file.path(folder_loc, "oc_post_sensitivity1.rds"))


# Sample Prior Sensitivity 1 ----------------------------------------------
oc_prior_sensitivity1 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                              data = model_objects_sensitivity1_priors_only,
                              init = init_func,
                              seed = 0,
                              chains = 4,
                              control = control_list)

write_rds(oc_prior_sensitivity1, file.path(folder_loc, "oc_prior_sensitivity1.rds"))

# Sample Posterior Sensitivity 2 ------------------------------------------
oc_post_sensitivity2 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                             data = model_objects_sensitivity2,
                             init = init_func,
                             seed = 0,
                             chains = 4,
                             control = control_list)

write_rds(oc_post_sensitivity2, file.path(folder_loc, "oc_post_sensitivity2.rds"))


# Sample Prior Sensitivity 2 ----------------------------------------------
oc_prior_sensitivity2 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                              data = model_objects_sensitivity2_priors_only,
                              init = init_func,
                              seed = 0,
                              chains = 4,
                              control = control_list)

write_rds(oc_prior_sensitivity2, file.path(folder_loc, "oc_prior_sensitivity2.rds"))
# Sample Posterior Sensitivity 3 ------------------------------------------
frac_carrs_prior3 <- mom_beta3(for_estimating$frac_carrs)

init_func3 <- function() list(frac_carrs = frac_carrs_prior3[1] / sum(frac_carrs_prior3),
                             frac_carrs_infec = frac_carrs_infec_prior[1] / sum(frac_carrs_infec_prior),
                             frac_infec_early = frac_infec_early_prior[1] / sum(frac_infec_early_prior))

oc_post_sensitivity3 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                             data = model_objects_sensitivity3,
                             init = init_func3,
                             seed = 0,
                             chains = 4,
                             control = control_list)

write_rds(oc_post_sensitivity3, file.path(folder_loc, "oc_post_sensitivity3.rds"))


# Sample Prior Sensitivity 3 ----------------------------------------------

oc_prior_sensitivity3 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                              data = model_objects_sensitivity3_priors_only,
                              init = init_func3,
                              seed = 0,
                              chains = 4,
                              control = control_list)

write_rds(oc_prior_sensitivity3, file.path(folder_loc, "oc_prior_sensitivity3.rds"))

# Sample Posterior Sensitivity 4 ------------------------------------------
frac_carrs_prior4 <- mom_beta4(for_estimating$frac_carrs)

init_func4 <- function() list(frac_carrs = frac_carrs_prior4[1] / sum(frac_carrs_prior4),
                              frac_carrs_infec = frac_carrs_infec_prior[1] / sum(frac_carrs_infec_prior),
                              frac_infec_early = frac_infec_early_prior[1] / sum(frac_infec_early_prior))

oc_post_sensitivity4 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                             data = model_objects_sensitivity4,
                             init = init_func4,
                             seed = 0,
                             chains = 4,
                             control = control_list)

write_rds(oc_post_sensitivity4, file.path(folder_loc, "oc_post_sensitivity4.rds"))


# Sample Prior Sensitivity 4 ----------------------------------------------

oc_prior_sensitivity4 <- stan(file = here("code", model_name, str_c(model_name, ".stan")),
                              data = model_objects_sensitivity4_priors_only,
                              init = init_func4,
                              seed = 0,
                              chains = 4,
                              control = control_list)

write_rds(oc_prior_sensitivity4, file.path(folder_loc, "oc_prior_sensitivity4.rds"))




}

#list.dirs(here("code", model_name, loc_name), recursive = F, full.names = F)

full_county_popsize <- 3175692
# Pick folder to build model from. Defaults to most recent end date with earliest start date
folder_name <- tibble(dirs = list.dirs(here("code", model_name, loc_name), recursive = F, full.names = F)) %>%
  filter(str_detect(dirs, "^\\d{4}-\\d{2}-\\d{2}_\\d{4}-\\d{2}-\\d{2}$")) %>% 
  separate(dirs, c("start", "end"), "_") %>%
  filter(end == max(end)) %>%
  #filter(start == "2020-04-06") %>%
  unite(col = dir, sep = "_") %>%
  pull(dir)

folder_loc <- here("code", model_name, loc_name, folder_name)

oc_post <- read_rds(file.path(folder_loc, "oc_post.rds"))
model_objects <- read_rds(file.path(folder_loc, "model_objects.rds"))

lumped_ochca_covid <- model_objects$lumped_ochca_covid
time_interval_in_days <- model_objects$time_interval_in_days

oc_prior <- read_rds(file.path(folder_loc, "oc_prior.rds"))

```

## Orange County, CA COVID-19 Situation Report, December 1st, 2020
#### Report period: `r str_c(format(model_objects$actual_first_day, "%b %d"), "-", format(model_objects$actual_last_day, "%b %d"), sep = " ")` (we don't use the most recent data due to reporting delays)

The goal of this report is to inform interested parties about dynamics of SARS-CoV-2 spread in Orange County, CA and to predict epidemic trajectories. 
Methodological details are provided below and in the accompanying [manuscript](https://arxiv.org/abs/2009.02654). We are also contributing to [COVID Trends by UC Irvine](https://www.stat.uci.edu/covid19/index.html) project that provides data visualizations of California County trends across time and space. 

```{r executive-summary, warning=FALSE, cache=TRUE}
R0_bci = oc_post %>% spread_draws(log_R0) %>% mutate(R0=exp(log_R0))%>% median_qi(R0) %>% mutate_if(is.numeric, signif, 2)

R_eff_bci <- spread_draws(oc_post, log_R0) %>%
  select(.draw, log_R0) %>%
  left_join(prep_epi_curves(oc_post, model_objects) %>%
              filter(g_text == "S",
                     t == model_objects$actual_last_day) %>%
              ungroup() %>%
              select(S = epi_curves, .draw)) %>%
  transmute(Reff = exp(log_R0) * S * model_objects$popsize / full_county_popsize) %>% 
  median_qi(Reff) %>% 
  mutate_if(is.numeric, signif, 2)

IFR_bci = oc_post %>% spread_draws(IFR) %>% median_qi(IFR) %>% mutate_if(is.numeric, signif, 2)

rho_death_bci = oc_post %>% spread_draws(rho_death) %>% median_qi(rho_death) %>% mutate_if(is.numeric, signif, 2)

ci_width <- 0.95
prepped_cumulative_inc <- prep_cumulative_inc(model_objects, oc_post)

death_inc_summaries <- prepped_cumulative_inc$posterior %>% 
  filter(usage == "train") %>%
  dplyr::select(-usage) %>% 
  group_by(t, name) %>% 
  median_qi(.width = ci_width) %>% 
  left_join(prepped_cumulative_inc$observed) %>% 
  mutate(ratio_l = .lower / observed) %>%
  mutate(ratio_u = .upper / observed) %>% 
  group_by(name) %>%
  filter(t == max(t)) %>%
  select(name, c(t, ratio_l, ratio_u, .lower, .upper, observed)) %>%
  mutate_if(is.numeric, signif, 2) %>%
  column_to_rownames(var = "name")
```


```{r process-calcat, warning=FALSE, cache=FALSE}
death_offset <- oc_data %>%
  filter(posted_date <= model_objects$actual_last_day) %>% 
  pull(new_deaths) %>% 
  sum()

data_for_calcat_re <- spread_draws(oc_post, log_R0) %>%
  select(.draw, log_R0) %>%
  left_join(prep_epi_curves(oc_post, model_objects) %>%
              filter(g_text == "S")) %>% 
  filter(t > model_objects$actual_last_day) %>% 
  select(date = t, log_R0, S = epi_curves) %>%
  group_by(date) %>% 
  transmute(Reff = exp(log_R0) * S * model_objects$popsize / full_county_popsize) %>% 
  median_qi() %>% 
  select(date, re_mean = Reff, re_CI95l = .lower, re_CI95u = .upper) %>% 
  right_join(tibble(date = seq(model_objects$actual_last_day, max(model_objects$dates_pp), 1))) %>% 
  pivot_longer(-date) %>% 
  arrange(date) %>%
  group_by(name) %>%
  mutate(id.value = na.approx(value, maxgap = 2)) %>% # linear interpolation for daily values
  pivot_wider(id_cols = date, names_from = name, values_from = id.value) %>% 
  select(date, re_mean, re_CI95l, re_CI95u)


data_for_calcat_mort <- oc_post %>%
  spread_draws(deaths_pp[i]) %>% 
  group_by(.draw) %>% 
  arrange(.draw, i) %>%
  mutate(date = model_objects$dates_pp[i]) %>% 
  filter(date > model_objects$actual_last_day) %>% 
  mutate(cummort = cumsum(deaths_pp) + death_offset) %>%
  ungroup() %>% 
  select(cummort, date) %>% 
  group_by(date) %>% 
  mean_qi() %>% 
  select(date, cummort_mean = cummort, cummort_CI95l = .lower, cummort_CI95u = .upper) %>% 
  right_join(tibble(date = seq(model_objects$actual_last_day, max(model_objects$dates_pp), 1))) %>% 
  pivot_longer(-date) %>% 
  arrange(date) %>%
  group_by(name) %>%
  mutate(id.value = na.approx(value, maxgap = 2)) %>% # linear interpolation for daily values
  pivot_wider(id_cols = date, names_from = name, values_from = id.value) %>% 
  select(date, cummort_mean, cummort_CI95l, cummort_CI95u)

data_for_calcat <- full_join(data_for_calcat_re, data_for_calcat_mort) %>% 
  drop_na() %>% 
  mutate(county = "Orange County", fips = 6059) %>% 
  select(date, county, fips, everything())

prev_data_for_calcat <- read_csv("data/data_for_calcat.csv")

write_csv(bind_rows(
  filter(prev_data_for_calcat, date %in% setdiff(prev_data_for_calcat$date, data_for_calcat$date)),
         data_for_calcat),
  path = "data/data_for_calcat.csv")
```

```{r latent-vs-observed-incidence, warning=FALSE, cache=TRUE}
ci_width <- c(0.5, 0.8, 0.95)
prepped_cumulative_inc <- map(prepped_cumulative_inc, ~mutate(., name = str_replace(name,
                                  "(Cumulative.+)",
                                  str_c("\\1\nSince ", format(model_objects$actual_first_day, "%b %d")))))
ggplot() +
  geom_lineribbon(data = prepped_cumulative_inc$posterior %>% 
                    group_by(t, name, usage) %>% 
                    median_qi(.width = ci_width),
  mapping = aes(t, value, ymin = .lower, ymax = .upper), size=1.5,
  color = "steelblue4") +
  geom_col(data = prepped_cumulative_inc$observed %>%
             filter(usage == "train"),
           mapping = aes(x = t, y = observed),
           fill = "black") +
  facet_wrap(. ~ name, scales = "free_y", strip.position="left") +
  ggtitle(str_c("OC latent & observed trajectories, posterior median &", str_c(percent(ci_width), collapse = ", "), "credible intervals", sep = " ")) +
  xlab("Date") +
  ylab(NULL) +
  theme(strip.background = element_blank(),
           strip.placement = "outside", 
        strip.text = element_text(size = 22),
        legend.position = c(0.092, 0.9), legend.title = element_text(size=15.5), legend.text = element_text(size=15.5), legend.background = element_rect(fill="transparent")) + 
  scale_fill_brewer(name="Credibility level", labels=str_c(percent(rev(ci_width))), guide = guide_legend(title.position = "top", direction = "horizontal")) +
  scale_x_date(breaks=c("21 day"), date_labels = "%b %d") +
  scale_y_continuous(labels = comma)
```

```{r ppc-forecasts, warning=FALSE, cache=TRUE}
ci_width <- c(0.5, 0.8, 0.95)
prepped_death <- prep_death(model_objects, oc_post, ci_width)

death_plot <- ggplot() +
  geom_lineribbon(data = prepped_death$posterior %>% 
                    group_by(t) %>% 
                    median_qi(.width = ci_width),
                  mapping = aes(t, deaths, ymin = .lower, ymax = .upper),
                  size = 1.5,
                  color = "steelblue4") +
  geom_point(data = prepped_death$observed,
             mapping = aes(t, observed), size=2,
             fill = "black") +
  #ggtitle("Deaths due to COVID-19") +
  xlab("Date") +
  ylab("Number of Deaths") +
  theme(legend.position = c(0.75, 0.8), legend.title = element_text(size=16), legend.text = element_text(size=16), legend.background = element_rect(fill="transparent")) + 
  scale_fill_brewer(name="Credibility level", labels=str_c(percent(rev(ci_width))), guide = guide_legend(title.position = "top", direction = "horizontal")) +
  scale_x_date(breaks=c("14 day"), date_labels = "%b %d")

prepped_pos <- prep_pos(model_objects, oc_post, ci_width)

test <- prepped_pos$posterior %>% 
  group_by(t) %>% 
  median_qi(.width = ci_width) 

positivity_plot <- ggplot() +
  geom_lineribbon(data = prepped_pos$posterior %>% 
                    group_by(t) %>% 
                    median_qi(.width = ci_width),
                  mapping = aes(t, percent_positive, ymin = .lower, ymax = .upper),
                  size = 1.5, color = "steelblue4") +
  geom_point(data = prepped_pos$observed,
             mapping = aes(t, observed), size=2,
             fill = "black") +
  #ggtitle("Fraction of Positive COVID-19 Tests") +
  xlab("Date") +
  ylab("Positivity Percent") +
  scale_y_continuous(labels = percent) +
  theme(legend.position = "none", legend.title = element_text(size=16), legend.text = element_text(size=16), legend.background = element_rect(fill="transparent")) + 
  scale_fill_brewer(name="Credibility level", labels=str_c(percent(rev(ci_width))), guide = guide_legend(title.position = "top", direction = "horizontal")) +
  scale_x_date(breaks=c("14 day"), date_labels = "%b %d")


death_pos <- death_plot + positivity_plot +   
  plot_annotation(title = 'Observed & Predicted Deaths and Positive Test Percent in 3 Day Periods', theme=theme(plot.title = element_text(hjust = 0.2)))

death_pos

```

```{r R0-IFR, echo=FALSE, cache=TRUE}
# If you want to plot R0 by istelf with custom lims, use this:
R0 <- oc_prior %>%
  spread_draws(log_R0) %>%
  mutate(R0=exp(log_R0))%>%
  dplyr::select(-starts_with(".")) %>%
  mutate(dist = "Prior") %>%
  bind_rows(oc_post %>%
              spread_draws(log_R0)%>% mutate(R0=exp(log_R0)) %>%
              dplyr::select(-starts_with(".")) %>%
              mutate(dist = "Posterior")) %>%
  ggplot(aes(x = R0, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  xlab(expression(R[0])) +
  coord_cartesian(xlim = c(0, 2.5))

IFR <- oc_prior %>%
  spread_draws(IFR) %>% 
  select(-starts_with(".")) %>% 
  mutate(dist = "Prior") %>% 
  bind_rows(oc_post %>%
              spread_draws(IFR) %>% 
              select(-starts_with(".")) %>% 
              mutate(dist = "Posterior"))  %>%
  ggplot(aes(x =IFR, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  #geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  xlab("Infection-to-fatality ratio (IFR)") +
  coord_cartesian(xlim = c(0, .03))

param_plot <- R0 + IFR +
    plot_annotation(title = expression(paste('Prior and Posterior Summaries for ', R[0], ' and IFR')),theme=theme(plot.title = element_text(hjust = 0.14)))

param_plot
```


```{r R0-over-time, cache=TRUE}
R0_over_time_plot <- tibble(file_path = fs::dir_ls(here("code", model_name, loc_name))) %>% 
  mutate(start_date = ymd(str_sub(file_path, start = -21, end = -12)),
    end_date = ymd(str_sub(file_path, start = -10))) %>% 
  mutate(fit_number = n():1) %>% 
  group_by(lubridate::month(end_date)) %>% 
  mutate(fit_number_month = row_number()) %>% 
  filter(fit_number %in% 1:4 | fit_number_month == 1) %>% 
  ungroup() %>% 
  select(file_path, start_date, end_date) %>% 
  filter(start_date > ymd("2020-06-16")) %>% # comment this out to include models not in archive
  # mutate(label = fct_inorder(glue("{format(start_date, format = '%b %e')}\n-\n{format(end_date, format = '%b %e')}"))) %>% 
  mutate(label = fct_inorder(glue("{format(start_date, format = '%b %e')} - {format(end_date, format = '%b %e')}"))) %>% 
  mutate(post = map(file_path, ~read_rds(path(., "oc_post.rds")))) %>% 
  mutate(tmp = map(post, ~spread_draws(., log_R0))) %>% 
  unnest(tmp) %>% 
  mutate(R0 = exp(log_R0)) %>%
  select(label, R0) %>% 
  ggplot(aes(label, R0)) +
  stat_eye(fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4") + 
  xlab("Date Range") +
  ylab(expression(R[0])) +
  geom_hline(yintercept = 1, linetype = "dashed") + 
plot_annotation(title = expression(paste('Historical estimates of the basic reproductive number, ', R[0], ', from previous reports')),theme=theme(plot.title = element_text(hjust = 0.4)))

R0_over_time_plot
```

## Summary (statements are made assuming 95% credibility levels)
- The number of reported cases (`r comma(death_inc_summaries["Cumulative Incidence","observed"])` in this period, shown as black bars in the top-middle plot above) underestimates the actual number of infections by a factor that ranges between `r death_inc_summaries["Cumulative Incidence","ratio_l"]` and `r death_inc_summaries["Cumulative Incidence","ratio_u"]`. 
This means that we estimate that the total number of infections which occurred between `r format(prepped_cumulative_inc$posterior[[1,"t"]]-2, '%B %d, %Y')` and `r format(death_inc_summaries["Cumulative Incidence","t"], '%B %d, %Y')` is between `r comma(death_inc_summaries["Cumulative Incidence",".lower"])` and `r comma(death_inc_summaries["Cumulative Incidence",".upper"])`.
- Prevalence (number of infectious individuals at any time point) is `r ifelse(median(filter(prepped_cumulative_inc$posterior, name == "Prevalence", t == min(t))$value) > median(filter(prepped_cumulative_inc$posterior, name == "Prevalence", t == max(t))$value), "decreasing", "increasing")` and projected to be inside the interval (`r prepped_cumulative_inc$posterior %>% filter(name == "Prevalence", t == max(t)) %>% pull(value) %>% quantile(c(0.025, 0.975)) %>% signif(digits=2) %>% scales::comma() %>% paste(collapse = ", ")`) on `r format(max(prepped_cumulative_inc$posterior$t), "%B %e")`.
- Deaths are also underreported, but not significantly so (black bars are a little below blue bands in the topleft plot above). Somewhere between `r signif(100*IFR_bci[".lower"],1)`% and `r signif(100*IFR_bci[".upper"],1)`% of all infections (not cases!) result in death.
- Basic reproductive number ($R_0$), defined as the average number of secondary infections one infectious individual produces in a completely susceptible population, is inside the interval (`r R0_bci[".lower"]`, `r R0_bci[".upper"]`). 
- Effective reproductive number ($R_e$), defined as its basic counterpart above, but allowing for some fraction of the population to be removed (recovered or deceased), as of `r format(death_inc_summaries["Cumulative Incidence","t"], '%B %d, %Y')` is inside the interval (`r R_eff_bci[".lower"]`, `r R_eff_bci[".upper"]`). **We want to keep $R_e < 1$ in order to control virus transmission**.

<hr style="height:3px;border-width:0;color:gray;background-color:gray">

## Abbreviated technical details (optional)

Our approach is based on fitting a mechanistic model of SARS-CoV-2 spread to multiple sources of surveillance data. 
A more fleshed out method description is in the [manuscript](https://arxiv.org/abs/2009.02654). 

### Model inputs
Our method takes three time series as input: daily new tests, case counts, and deaths. However, we find daily resolution to be too noisy due to delay in testing reports, weekend effect, etc. So we aggregated/binned the three types of counts in 3 day intervals. These aggregated time series are shown below.

```{r model-inputs, cache=TRUE}
lumped_ochca_covid %>% 
  dplyr::select(-start_date) %>% 
  pivot_longer(-end_date) %>% 
  mutate(name = str_remove(name, "new_") %>% str_to_title()) %>% 
  mutate(name = fct_relevel(name, c("Tests", "Cases", "Deaths"))) %>% 
  ggplot(aes(end_date, value)) +
  geom_line(size=1) +
  geom_point(size=3) +
  facet_wrap(. ~ name, scales = "free_y") +
  xlab("Date") +
  ylab("Count") +
  scale_y_continuous(labels = comma) +
  ggtitle("Orange County, CA data", subtitle = str_c("Counts binned into", as.integer(lumped_ochca_covid$end_date[2] - lumped_ochca_covid$end_date[1]), "day periods", sep = " ")) + scale_x_date(breaks=c("10 day"), date_labels = "%b %d", )+
  theme(strip.text = element_text(size = 20))
```

### Model structure
We assume that all individuals in Orange County, CA can be split into 6 compartments: S = susceptible individuals, E = infected, but not yet infectious individuals, $\text{I}_\text{e}$ = individuals at early stages of infection, $\text{I}_\text{p}$ = individuals at progressed stages of infection (assumed 20% less infectious than individuals at the early infection stage), R = recovered individuals, D = individuals who died due to COVID-19. Possible progressions of an individual through the above compartments are depicted in the diagram below.  


```{r Model Structure, out.width = "60%"}
knitr::include_graphics("assets/model_figure.svg", error=FALSE)
```

Mathematically, we assume that dynamics of the proportions of individuals in each compartment follow a set of ordinary differential equations corresponding to the above diagram. These equations are controlled by the following parameters: 

- Basic reproductive number ($R_0$)
- mean duration of the latent period
- mean duration of the early infection period
- mean duration of the progressed infection period
- probability of transitioning from progressed infection to death, rather than to recovery (i.e., IFR)

We fit this model to data by assuming that case counts are noisy realizations of the actual number of individuals progressing from $\text{I}_\text{e}$ compartment to $\text{I}_\text{p}$ compartment. 
Similarly we assume that observed deaths are noisy realizations of the actual number of individuals progressing from $\text{I}_\text{p}$ compartment to $\text{D}$ compartment. 
*A priori*, we assume that death counts are significantly less noisy than case counts.
We use a Bayesian estimation framework, which means that all estimated quantities receive credible intervals (e.g., 80% or 95% credible intervals). 
Width of these credible intervals encode the amount of uncertainty that we have in the estimated quantities. 





<!-- # Planned model extensions  -->

<!-- * We plan to add H = hospitalized and ICU compartments to the model. This extension will allow us to inform the model with hospitalization and occupied ICU beds time series, in addition to already used test, case, and death counts. Equally importantly, we will be able to forecast hospital and ICU occupancy.  -->
<!-- * We will allow $R_0$ to vary with time and will attempt to estimate its changes non-parameterically. -->
<!-- * We would like to add age structure into the model. This is not a high priority at the moment, because we suspect we don't have enough data to estimate additional parameters that will need to be added to the model. -->

<hr style="height:3px;border-width:0;color:gray;background-color:gray">

## Appendix

### Sensitivity to Prior for $R_0$

We examine how sensitive our conclusions about $R_0$ are to our prior assumptions by repeating estimation of all model parameters under different priors for this parameter. 
The priors are listed in the titles of the figures.
Although the prior distribution of $R_0$ does have some effect on its posterior (as it should), our results and conclusions are not too sensitive to a particular specification of this prior. 


```{r reading-sensitivity-objects, echo=FALSE, cache=TRUE}
#Read in
#sensitivity1 is lognorm(log(2.5), .5)
#sensitivity2 is lognorm(log(.53), .78)
#folder_loc <- "Models/simple_mod_for_simulation/2020-03-15_2020-04-28"
oc_post_sensitivity1<- read_rds(file.path(folder_loc, "oc_post_sensitivity1.rds"))
oc_prior_sensitivity1 <- read_rds(file.path(folder_loc, "oc_prior_sensitivity1.rds"))

oc_post_sensitivity2<- read_rds(file.path(folder_loc, "oc_post_sensitivity2.rds"))
oc_prior_sensitivity2 <- read_rds(file.path(folder_loc, "oc_prior_sensitivity2.rds"))

oc_post_sensitivity3 <- read_rds(file.path(folder_loc, "oc_post_sensitivity3.rds"))
oc_prior_sensitivity3 <- read_rds(file.path(folder_loc, "oc_prior_sensitivity3.rds"))

oc_post_sensitivity4<- read_rds(file.path(folder_loc, "oc_post_sensitivity4.rds"))
oc_prior_sensitivity4 <- read_rds(file.path(folder_loc, "oc_prior_sensitivity4.rds"))


prepped_cumulative_inc_s1 <- prep_cumulative_inc(model_objects, oc_post_sensitivity1)
prepped_cumulative_inc_s2 <- prep_cumulative_inc(model_objects, oc_post_sensitivity2)
prepped_cumulative_inc_s3 <- prep_cumulative_inc(model_objects, oc_post_sensitivity3)
prepped_cumulative_inc_s4 <- prep_cumulative_inc(model_objects, oc_post_sensitivity4)
```

```{r, echo=FALSE, cache=TRUE}
ci_width <- c(0.5, 0.8, 0.95)

ggplot() +
  geom_lineribbon(data = prepped_cumulative_inc_s1$posterior %>%
  group_by(t, name, usage) %>%
  median_qi(.width = ci_width),
  mapping = aes(t, value, ymin = .lower, ymax = .upper),
  color = "steelblue4") +
  geom_col(data = prepped_cumulative_inc_s1$observed %>% 
                    filter(usage == "train"),
           mapping = aes(x = t, y = observed),
           fill = "black") +
  facet_wrap(. ~ name, scales = "free_y", strip.position = "left") +
  ggtitle("Latent & observed cumulative events R0~lognorm(log(2.5), 0.5)",
          subtitle = str_c("Posterior Median &", str_c(percent(ci_width), collapse = ", "), "credible intervals", sep = " ")) +
  xlab("Date") +
  ylab(NULL) +
  theme(strip.background = element_blank(),
           strip.placement = "outside", strip.text = element_text(size = 22),
        legend.position = c(0.092, 0.9), legend.title = element_text(size=15.5), legend.text = element_text(size=15.5), legend.background = element_rect(fill="transparent")) +
  scale_fill_brewer(name="Credibility level", labels=str_c(percent(rev(ci_width))), guide = guide_legend(title.position = "top", direction = "horizontal")) +
  scale_x_date(breaks=c("21 day"), date_labels = "%b %d") +
  scale_y_continuous(labels = comma)
```

```{r, warnings=FALSE, cache=TRUE}
## If you want to plot R0 by istelf with custom lims, use this:
R0_s1 <- oc_prior_sensitivity1 %>%
  spread_draws(log_R0) %>%
  mutate(R0=exp(log_R0))%>%
  select(-starts_with(".")) %>%
  mutate(dist = "Prior") %>%
  bind_rows(oc_post_sensitivity1 %>%
              spread_draws(log_R0) %>%
              mutate(R0=exp(log_R0))%>%
              select(-starts_with(".")) %>%
              mutate(dist = "Posterior")) %>%
  ggplot(aes(x = R0, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  xlab(expression(R[0])) +
  coord_cartesian(xlim = c(0, 2.5))
# and substitute


IFR_s1 <- oc_prior_sensitivity1%>%
  spread_draws(IFR) %>% 
  select(-starts_with(".")) %>% 
  mutate(dist = "Prior") %>% 
  bind_rows(oc_post_sensitivity1 %>%
              spread_draws(IFR) %>% 
              select(-starts_with(".")) %>% 
              mutate(dist = "Posterior"))  %>%
  ggplot(aes(x =IFR, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  #geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  coord_cartesian(xlim = c(0, .05))

param_plot_s1 <- R0_s1 + IFR_s1 +
    plot_annotation(title = 'Prior and Posterior summaries for R0 and IFR',
                    subtitle= 'R0 ~ lognormal(log(2.5), .5)')

param_plot_s1

```

```{r, warning=FALSE, cache=TRUE}
ci_width <- c(0.5, 0.8, 0.95)

ggplot() +
  geom_lineribbon(data = prepped_cumulative_inc_s2$posterior %>%
                    group_by(t, name, usage) %>%
                    median_qi(.width = ci_width),
                  mapping = aes(t, value, ymin = .lower, ymax = .upper),
                  color = "steelblue4") +
  geom_col(data = prepped_cumulative_inc_s2$observed %>% 
                    filter(usage == "train"),
           mapping = aes(x = t, y = observed),
           fill = "black") +
  facet_wrap(. ~ name, scales = "free_y", strip.position = "left") +
  ggtitle("Latent & observed cumulative events R0~lognorm(log(0.53), 0.78)",
          subtitle = str_c("Posterior Median &", str_c(percent(ci_width), collapse = ", "), "credible intervals", sep = " ")) +
  xlab("Date") +
  ylab(NULL) +
  theme(strip.background = element_blank(),
           strip.placement = "outside", strip.text = element_text(size = 22),
        legend.position = c(0.092, 0.9), legend.title = element_text(size=15.5), legend.text = element_text(size=15.5), legend.background = element_rect(fill="transparent")) +
  scale_fill_brewer(name="Credibility level", labels=str_c(percent(rev(ci_width))), guide = guide_legend(title.position = "top", direction = "horizontal")) +
  scale_x_date(breaks=c("21 day"), date_labels = "%b %d") +
  scale_y_continuous(labels = comma)

```

```{r, echo=FALSE, cache=TRUE}
R0_s2 <- oc_prior_sensitivity2 %>%
  spread_draws(log_R0) %>%
  mutate(R0=exp(log_R0))%>%
  select(-starts_with(".")) %>%
  mutate(dist = "Prior") %>%
  bind_rows(oc_post_sensitivity2 %>%
              spread_draws(log_R0) %>%
              mutate(R0=exp(log_R0))%>%
              select(-starts_with(".")) %>%
              mutate(dist = "Posterior")) %>%
  ggplot(aes(x = R0, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  xlab(expression(R[0])) +
  coord_cartesian(xlim = c(0, 2.5))
# and substitute


IFR_s2 <- oc_prior_sensitivity2%>%
  spread_draws(IFR) %>% 
  select(-starts_with(".")) %>% 
  mutate(dist = "Prior") %>% 
  bind_rows(oc_post_sensitivity2 %>%
              spread_draws(IFR) %>% 
              select(-starts_with(".")) %>% 
              mutate(dist = "Posterior"))  %>%
  ggplot(aes(x =IFR, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  #geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  coord_cartesian(xlim = c(0, .05))

param_plot_s2 <- R0_s2 + IFR_s2 +
    plot_annotation(title = expression(paste('Prior and Posterior summaries for ', R[0], ' and IFR')),
                    subtitle= 'R0 ~ lognormal(log(.53), .78)')

param_plot_s2
```


### Sensitivity to prior for fraction initially infected
We examine how sensitive our conclusions about $R_0$ are to our prior assumptions by repeating estimation of all model parameters under different priors for the parameter controlling how many people are infected initially. This prior changes depending on the time period, so we adjust by changing the prior mean to be twice as large or one half as large as the default prior. 

<!-- sensitivity 3 -->

```{r, echo=FALSE, cache=TRUE}
ci_width <- c(0.5, 0.8, 0.95)

ggplot() +
  geom_lineribbon(data = prepped_cumulative_inc_s3$posterior %>%
  group_by(t, name, usage) %>%
  median_qi(.width = ci_width),
  mapping = aes(t, value, ymin = .lower, ymax = .upper),
  color = "steelblue4") +
  geom_col(data = prepped_cumulative_inc_s3$observed %>% 
                    filter(usage == "train"),
           mapping = aes(x = t, y = observed),
           fill = "black") +
  facet_wrap(. ~ name, scales = "free_y", strip.position = "left") +
  ggtitle("Latent & observed cumulative events: Mean Initial Infected Twice as Large",
          subtitle = str_c("Posterior Median &", str_c(percent(ci_width), collapse = ", "), "credible intervals", sep = " ")) +
  xlab("Date") +
  ylab(NULL) +
  theme(strip.background = element_blank(),
           strip.placement = "outside", strip.text = element_text(size = 22),
        legend.position = c(0.092, 0.9), legend.title = element_text(size=15.5), legend.text = element_text(size=15.5), legend.background = element_rect(fill="transparent")) +
  scale_fill_brewer(name="Credibility level", labels=str_c(percent(rev(ci_width))), guide = guide_legend(title.position = "top", direction = "horizontal")) +
  scale_x_date(breaks=c("21 day"), date_labels = "%b %d") +
  scale_y_continuous(labels = comma)
```

```{r, warnings=FALSE, cache=TRUE}
R0_s3 <- oc_prior_sensitivity3 %>%
  spread_draws(log_R0) %>%
  mutate(R0=exp(log_R0))%>%
  select(-starts_with(".")) %>%
  mutate(dist = "Prior") %>%
  bind_rows(oc_post_sensitivity3 %>%
              spread_draws(log_R0) %>%
              mutate(R0=exp(log_R0))%>%
              select(-starts_with(".")) %>%
              mutate(dist = "Posterior")) %>%
  ggplot(aes(x = R0, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  xlab(expression(R[0])) +
  coord_cartesian(xlim = c(0, 2.5))
# and substitute


IFR_s3 <- oc_prior_sensitivity3%>%
  spread_draws(IFR) %>% 
  select(-starts_with(".")) %>% 
  mutate(dist = "Prior") %>% 
  bind_rows(oc_post_sensitivity3 %>%
              spread_draws(IFR) %>% 
              select(-starts_with(".")) %>% 
              mutate(dist = "Posterior"))  %>%
  ggplot(aes(x =IFR, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  #geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  coord_cartesian(xlim = c(0, .05))

param_plot_s3 <- R0_s3 + IFR_s3 +
    plot_annotation(title = expression(paste('Prior and Posterior summaries for ', R[0], ' and IFR')),
                    subtitle= 'Mean Initial Infected Twice as Large')

param_plot_s3
```

<!-- sensitivity 4 -->

```{r, echo=FALSE, cache=TRUE}
ci_width <- c(0.5, 0.8, 0.95)

ggplot() +
  geom_lineribbon(data = prepped_cumulative_inc_s4$posterior %>%
  group_by(t, name, usage) %>%
  median_qi(.width = ci_width),
  mapping = aes(t, value, ymin = .lower, ymax = .upper),
  color = "steelblue4") +
  geom_col(data = prepped_cumulative_inc_s4$observed %>% 
                    filter(usage == "train"),
           mapping = aes(x = t, y = observed),
           fill = "black") +
  facet_wrap(. ~ name, scales = "free_y", strip.position = "left") +
  ggtitle("Latent & observed cumulative events: Mean Initial Infected Half as Large",
          subtitle = str_c("Posterior Median &", str_c(percent(ci_width), collapse = ", "), "credible intervals", sep = " ")) +
  xlab("Date") +
  ylab(NULL) +
  theme(strip.background = element_blank(),
           strip.placement = "outside", strip.text = element_text(size = 22),
        legend.position = c(0.092, 0.9), legend.title = element_text(size=15.5), legend.text = element_text(size=15.5), legend.background = element_rect(fill="transparent")) +
  scale_fill_brewer(name="Credibility level", labels=str_c(percent(rev(ci_width))), guide = guide_legend(title.position = "top", direction = "horizontal")) +
  scale_x_date(breaks=c("21 day"), date_labels = "%b %d") +
  scale_y_continuous(labels = comma)
```

```{r, warnings=FALSE, cache=TRUE}
R0_s4 <- oc_prior_sensitivity4 %>%
  spread_draws(log_R0) %>%
  mutate(R0=exp(log_R0))%>%
  select(-starts_with(".")) %>%
  mutate(dist = "Prior") %>%
  bind_rows(oc_post_sensitivity4 %>%
              spread_draws(log_R0) %>%
              mutate(R0=exp(log_R0))%>%
              select(-starts_with(".")) %>%
              mutate(dist = "Posterior")) %>%
  ggplot(aes(x = R0, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  xlab(expression(R[0])) +
  coord_cartesian(xlim = c(0, 2.5))
# and substitute


IFR_s4 <- oc_prior_sensitivity4%>%
  spread_draws(IFR) %>% 
  select(-starts_with(".")) %>% 
  mutate(dist = "Prior") %>% 
  bind_rows(oc_post_sensitivity4 %>%
              spread_draws(IFR) %>% 
              select(-starts_with(".")) %>% 
              mutate(dist = "Posterior"))  %>%
  ggplot(aes(x =IFR, y = dist), fill = stat(x > 1.0)) +
  stat_halfeye(normalize= "xy", fill="lightskyblue1", slab_color ="lightskyblue4", color="dodgerblue4", slab_size = 0.5) +
  #geom_vline(xintercept = 1.0, linetype = "dashed") +
  ylab(NULL) +
  coord_cartesian(xlim = c(0, .05))

param_plot_s4 <- R0_s4 + IFR_s4 +
    plot_annotation(title = expression(paste('Prior and Posterior summaries for ', R[0], ' and IFR')),
                    subtitle= 'Mean Initial Infected Half as Large')

param_plot_s4
```

Last updated on `r Sys.Date()`.