suppressPackageStartupMessages({
library(tidyverse)
library(ggpubr)
library(lubridate)
library(knitr)
library(data.table)
source("scripts/data_preproc/preproc_functions.R")
source("scripts/mixture_analysis/functions.R")
})
options(digits=4)
knitr::opts_chunk$set(cache=TRUE)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
loaded_data <- new.env()
load("results/region-sparse-oneintervention-mixture-gamma-prior-latest.Rdata",
v=T, envir=loaded_data)
## Loading objects:
## fit
## out
## grouped_dataset
## covariates
## stan_data
stanfit$alpha %>%
as.data.frame() %>%
set_names(paste("alpha", c(1,2))) %>%
pivot_longer(everything(), names_to="alpha", values_to="value") %>%
ggplot() +
geom_histogram(aes(value, fill=alpha), position="dodge")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 2.5% 25% 50% 75% 97.5%
## [1,] 1.458 1.533 1.566 1.596 1.651
## [2,] 1.667 1.744 1.792 1.842 1.941
total_pop <-
read.csv("data/region_wise_IFR.csv", sep=";") %>%
prepare_regions_pop_ifr() %>%
pull(pop) %>% sum()
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(c(0.025, 0.25, 0.5, 0.75, 0.975))
df %>% filter(date == LOCKDOWN_END) %>%
mutate(total_cases_prop=total_cases / total_pop) %>%
group_by(sample_id) %>%
summarise_at(vars(total_cases, total_cases_prop, total_deaths, new_cases), sum) %>%
pivot_longer(-sample_id, names_to="variable", values_to="value") %>%
group_by(variable) %>%
summarise_at(vars(value), q_funs)
df %>% filter(date %in% c(LOCKDOWN, LOCKDOWN - 1)) %>%
group_by(sample_id, date) %>%
summarise_at(vars(Rt), mean) %>%
group_by(date) %>%
summarise_at(vars(Rt), q_funs)
(Rt_foldchanges <- df %>%
mutate(Rt_foldchange = exp(-alpha)) %>%
group_by(region, region_id, sample_id) %>%
summarise(
Rt_foldchange = unique(Rt_foldchange),
R0 = unique(R0)
) %>%
group_by(region, region_id) %>%
summarise_at(vars(Rt_foldchange,R0), median)
)
## `summarise()` regrouping output by 'region', 'region_id' (override with `.groups` argument)
(
Rt_fc_plt <- df %>%
mutate(Rt_foldchange = exp(-alpha)) %>%
group_by(region, region_id, sample_id) %>%
summarise(Rt_foldchange = unique(Rt_foldchange)) %>%
ggplot(aes(fct_rev(region))) +
geom_boxplot(aes(y=Rt_foldchange)) +
scale_x_discrete("") +
scale_y_continuous("Rt fold-change") +
coord_flip() +
theme_pubr() +
theme(plot.margin = unit(c(0.5,1,0.5,0), "cm"))
)
## `summarise()` regrouping output by 'region', 'region_id' (override with `.groups` argument)
df %>%
mutate(squared_err = (deaths_pred - deaths)^2) %>%
group_by(sample_id) %>%
summarise_at(vars(squared_err), sum, na.rm=T) %>%
pull(squared_err) %>%
quantile(c(0.025, 0.25, 0.5, 0.75, 0.975))
## 2.5% 25% 50% 75% 97.5%
## 193776 233693 257950 283389 345351
df %>%
group_by(sample_id, region, region_id) %>%
mutate(
deaths= cumsum(replace_na(deaths, 0)),
predicted_deaths_csum = cumsum(deaths_pred)) %>%
group_by(date, sample_id) %>%
summarise_at(vars(deaths, starts_with("predicted")), sum, na.rm=T) %>%
group_by(date, deaths) %>%
summarise(
predicted_lower= quantile(predicted_deaths_csum, 0.025),
predicted_median= quantile(predicted_deaths_csum, 0.5),
predicted_upper= quantile(predicted_deaths_csum, 0.975)
) %>%
ggplot(aes(x=date)) +
geom_line(aes(y=deaths)) +
geom_line(aes(y=predicted_median), col="red", lty="dashed") +
geom_ribbon(aes(ymin=predicted_lower, ymax=predicted_upper), alpha=0.5)
## `summarise()` regrouping output by 'date' (override with `.groups` argument)
df %>%
filter(region_id == 8) %>%
group_by(date, deaths) %>%
summarise(
pred_lower = quantile(deaths_pred, 0.025),
pred_med = quantile(deaths_pred, 0.5),
pred_upper = quantile(deaths_pred, 0.975)) %>%
ggplot(aes(x=date)) +
geom_col(aes(y=deaths)) +
geom_line(aes(y=pred_med)) +
geom_ribbon(aes(ymin=pred_lower, ymax=pred_upper), fill="deepskyblue", alpha=0.3)
## `summarise()` regrouping output by 'date' (override with `.groups` argument)
## Warning: Removed 22 rows containing missing values (position_stack).