suppressPackageStartupMessages({
library(tidyverse)
library(lubridate)
library(pbapply)
library(data.table)
library(knitr)
knitr::opts_chunk$set(cache=TRUE)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
source("scripts/mixture_analysis/functions.R")
})
options(digits=3)
pboptions(type="timer")
LOCKDOWN <- ymd("2020-03-17")
LOCKDOWN_END <- ymd("2020-05-11")
loaded_data <- new.env()
load("results/region-sparse-oneintervention-mixture-gamma-prior-latest.Rdata",
envir=loaded_data)
dataset <- loaded_data$grouped_dataset %>% select(-t)
stanfit <- loaded_data$out
On commence par calculer la log posterior associée à chaque observation (par échantillon/ condition / jour / region). Le résultat est filtré seulement pour réduire la taille du doc.
(df <- extract_mixture_dataset(dataset, stanfit, conditions=c(1,2)) %>%
filter(region_id <=2, sample_id <= 20))
Calcul d’une posterior par échantillon / condition / region en intégrant sur l’ensemble des jours observés, et en pondérant par la valeur de \(\theta\) échantillonnée. La colonne posterior
correspondrait à \(\theta^{condition}_{échantillon} f(R^{condition}_{region, échantillon})\) ; elle est exprimée comme proba à partir de là.
(posteriors <- df %>%
group_by(sample_id, region, region_id, condition) %>%
summarise(posterior=exp(sum(log(theta) + posterior, na.rm=T))))
## `summarise()` regrouping output by 'sample_id', 'region', 'region_id' (override with `.groups` argument)
On assigne la posterior “globale” par échantillon/condition/région à l’ensemble des observations
(mixture_posteriors <- df %>%
select(-posterior) %>%
inner_join(posteriors, by=c("region", "region_id", "condition", "sample_id")))
Pondération des estimations entre conditions par leur posterior. C’est un calcul par groupe, un groupe= 1 échantillon, jour, région. Un groupe contient une paire d’estimations pour chaque variable qui nous intéresse (comme deaths_pred
), associée à une paire de posterior. Chaque paire est pondérée par les posteriors des deux conditions.
# Using data.table for speed
(mixture_decoded <-
as.data.table(mixture_posteriors)[ ,
lapply(.SD, function(value){sum(value * posterior) / sum(posterior)}),
by=.(region, region_id, day_id, date, sample_id, deaths),
.SDcols=c(
"deaths_mean", "deaths_pred", "Rt", "alpha", "new_cases",
"total_cases", "total_deaths")
])