suppressPackageStartupMessages({
  library(tidyverse)
  library(ggpubr)
  library(lubridate)
  library(knitr)
  knitr::opts_chunk$set(cache=TRUE)
  knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
  source("scripts/plots/utils/utils.R")
  source("scripts/data_preproc/preproc_functions.R")
})
options(digits=3)
LOCKDOWN <- ymd("2020-03-17")
Rt_extractor <- function(region_data, keys){
  lockdown_day_id <-  LOCKDOWN - min(region_data$date) + 1
  
  stanfit$Rt_adj[,c(1, lockdown_day_id),keys$region_id] %>%
    as.data.frame() %>%
    set_names(c("R0", "Rt lockdown")) %>%
    pivot_longer(cols=everything(), names_to = "variable", values_to = "value")
}

1 Reproductive number

1.1 Base model vs week-ends model

data_WE <- new.env()
load("results/region-sparse-weekends-latest.Rdata", verbose=T, envir = data_WE)
## Loading objects:
##   fit
##   out
##   grouped_dataset
##   covariates
##   stan_data
stanfit <- data_WE$out

Rt_WE_df <- data_WE$grouped_dataset %>% 
  group_modify(Rt_extractor) %>%
  mutate(model="Week-ends")


data_m0 <- new.env()
load("results/region-sparse-lockdown-latest.Rdata", verbose=T, envir = data_m0)
## Loading objects:
##   fit
##   out
##   grouped_dataset
##   covariates
##   stan_data
stanfit <- data_m0$out

Rt_m0_df <- data_m0$grouped_dataset %>% 
  group_modify(Rt_extractor) %>%
  mutate(model='Base')

Rt_df <- rbind(Rt_WE_df, Rt_m0_df) %>% arrange(region, variable) %>%
  mutate(
    model = factor(model),
    variable= factor(variable)
  )

(
  plt <- Rt_df %>% 
    ggplot(aes(factor(region) %>% fct_rev())) +
    geom_boxplot(aes(y=value, group=interaction(region, model),fill=model), outlier.size = 0.5, outlier.alpha = 0.5) +
    facet_wrap(~variable, scales = "free_x")+
    coord_flip() +
    ylab('Estimated reproductive number') +
    theme_pubr() +
    theme(axis.title.y=element_blank())
)

cowplot::save_plot('report/SupplementaryData/figures/comparison_WE_base.png', plt, scale=1.2)