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")
    ])