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")
}
Reproductive number
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)