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)
invisible(Sys.setlocale("LC_ALL", locale="en_US.utf8"))
LOCKDOWN <- ymd("2020-03-17")
LOCKDOWN_END <- ymd("2020-05-11")
suppressMessages(
grouped_dataset <-
read_csv("results/merge_spf_opencovid_agences_regionales_sante/mergedDataReadyToAnalyze.csv") %>%
prepare_french_region_dataset() %>%
group_by(region, region_id) %>%
group_modify(function (region_data, keys) {
message(keys$region %>% str_to_upper())
message(" + PADDING : ")
region_data <- pad_region_data(region_data)
message(" + SLICING : ")
region_data <- slice_from_onset(region_data)
})
)
obs_total_deaths <- sum(grouped_dataset$deaths)
mixture <- data.table::fread(
file="results/simulated_datasets/mixture_lockdown_shift_decoded_posterior.csv",
colClasses = c(date="Date", lockdown_day="Date")
) %>%
filter(date <= LOCKDOWN_END) %>%
mutate(
model="mixture",
sim_id = dense_rank(sample_id)
) %>%
select(
region, region_id, sim_id, sample_id, day_id, date,
deaths_mean, deaths_pred, Rt, new_cases, lockdown_day, model
) %>%
as_tibble()
base_results <- new.env()
load("results/simulated_datasets/lockdown_shift.Rdata", env=base_results)
simulated_dataset <- base_results$simulated_dataset %>%
ungroup() %>%
filter(date <= LOCKDOWN_END) %>%
mutate(model="base") %>%
select(
region, region_id, sim_id, sample_id, day_id, date,
deaths_mean, deaths_pred, Rt, new_cases, lockdown_day, model
)
merged_dataset <- rbind(simulated_dataset, mixture) %>%
as.data.frame() %>%
filter(
date <= LOCKDOWN_END,
lockdown_day %>% between(LOCKDOWN-7, LOCKDOWN + 7)
) %>%
group_by(lockdown_day, sim_id, model) %>%
summarise(total_deaths=sum(na.omit(deaths_pred))) %>%
mutate(actual = lockdown_day == LOCKDOWN) %>%
ungroup()
## `summarise()` regrouping output by 'lockdown_day', 'sim_id' (override with `.groups` argument)
(
plt_shift <- merged_dataset %>%
ggplot(aes(x=lockdown_day)) +
geom_hline(
aes(yintercept=obs_total_deaths, lty="observed"),
col="orangered", lwd=0.4) +
geom_boxplot(
aes( y=total_deaths, group=interaction(lockdown_day, model), fill=model),
lwd=0.3) +
scale_x_date("Lockdown start", breaks="2 days", date_labels="%b %d") +
scale_y_continuous("Total mortality") +
scale_fill_discrete("Model") +
scale_linetype_manual(NULL, values=c("dashed")) +
theme_pubr()
)
deaths_lockdown_diff <- function(dataset, obs_total_deaths){
dataset %>%
group_by(lockdown_day, sample_id) %>%
summarise(total_deaths = sum(deaths_pred, na.rm=T)) %>%
group_by(lockdown_day) %>%
summarise_at(vars(total_deaths), median) %>%
# total_deaths_factual <- d %>% filter(lockdown_day == LOCKDOWN) %>% pull(total_deaths)
mutate(
diff = total_deaths - obs_total_deaths,
diff_perc = diff / obs_total_deaths)
}
## `summarise()` regrouping output by 'lockdown_day' (override with `.groups` argument)
## `summarise()` regrouping output by 'lockdown_day' (override with `.groups` argument)
q <- c(0.025, 0.25, 0.5, 0.75, 0.975)
q_funs <- map(q, ~partial(quantile, probs=.x, na.rm=TRUE)) %>% set_names(paste0(q*100, "%"))
(
total_deaths_1303 <- merged_dataset %>% filter(lockdown_day == ymd("2020-03-13")) %>%
group_by(model) %>%
summarise_at(vars(total_deaths), q_funs)
)