S6.Rmd

---
title: "S6"
output: 
  html_document:
    toc: true
    toc_float:
      smooth_scroll: FALSE
    number_sections: true
---

This appendix aims to illustrate the performance of the alternative 
parameterisation under various transmissibility levels ($\Re_0$) and mean 
generation times. For simplicity, we restrict the $SE^iI^jR$ to four instances
($i = 1$, and $j = \{1,2,3,4\}$) for data generation and fitting.


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)

library(cmdstanr)
library(dplyr)
library(Metrics)
library(purrr)
library(readr)
library(readsdr)
library(stringr)
library(tictoc)
library(tidyr)

source("./R_scripts/Inference.R")
source("./R_scripts/formulas.R")
source("./R_scripts/model_selection.R")
source("./R_scripts/plots.R")
source("./R_scripts/posterior_utils.R")
source("./R_scripts/synthetic_data.R")
source("./R_scripts/tables.R")
source("./R_scripts/utils.R")
```

# Overview

This sensitivity analysis consists of testing the alternative parameterisation
under **four** additional scenarios. The first scenario corresponds to that 
presented in **S3**. By *scenario*, we refer to the particular 
configuration of $\Re_0$ and the mean generation time ($\tau$) in the $SE^iI^jR$
framework to produce synthetic data. Since, $\tau$ varies as the distribution of
the infectious period ($j$) changes, we identify each scenario by the mean 
generation time obtained from an exponentially-distributed infectious period 
(denoted by $\tau_e$). For each instance and data-fidelity level, we generate
20 synthetic incidence reports. Recall that fidelity levels correspond to two
configurations of the Negative Binomial measurement model: high fidelity
($\phi^{-1} = 0$) & low fidelity ($\phi^{-1} = \frac{1}{3}$). As a result, we 
produce *160* synthetic incidence reports per scenario. Then, we fit four 
candidates per incidence report, where we assume the appropriate 
measurement model.

```{r}
scenario_df <- data.frame(Scenario = 1:5,
                          R0       = c(2.5, 2.5, 2.5, 9, 17),
                          tau      = c(4, 8, 13, 4, 4))

table_scenarios(scenario_df)
```


```{r}
params_df <- read_csv("./data/param_values.csv", show_col_types = FALSE)

# Population size
N_val         <- params_df |> filter(name == "N") |> pull(value)
# Reporting rate
rho_val       <- params_df |> filter(name == "rho") |> pull(value)

n_ds        <- 20
m_val       <- 1
n_val       <- 1:4
inv_phi_val <- 1/3

i_val   <- 1
j_val   <- 1:4
M_j     <- j_val

inv_phi_pois <- 0
inv_phi_ovds <- 1/3 # Overdispersed
```

```{r inferece_params}
n_chains          <- 4
samples_per_chain <- 1000

M_i        <- 1
M_j        <- 1:4
n_params   <- 5

all_nbin_files  <- readRDS("./Stan_files/Inference/Three_params_alt/nbin/meta_info.rds")
fltr            <- map_lgl(all_nbin_files, \(fl_list) fl_list$E_ord == M_i)
info_files_nbin <- all_nbin_files [fltr]

all_pois_files  <- readRDS("./Stan_files/Inference/Three_params_alt/pois/meta_info.rds")
fltr            <- map_lgl(all_pois_files, \(fl_list) fl_list$E_ord == M_i)
info_files_pois <- all_pois_files [fltr]


files_trad_pois <- readRDS("./Stan_files/Inference/Three_params/pois/meta_info.rds")
fltr            <- map_lgl(files_trad_pois, \(fl_list) fl_list$E_ord == M_i)
info_trad_pois  <- files_trad_pois[fltr]
```

# Scenario 2

In comparison to Scenario 1, we increase the reference mean generation time 
($\tau_e$) from **4** to **8**.

## Synthetic data

```{r}
inv_sigma_val <- 3
inv_gamma_val <- 5 # Medium
beta_val      <- 0.5 # R0 2.5
R0_val        <- beta_val * inv_gamma_val
ts_val        <- 100
```


```{r}
set.seed(1600)

y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_ovds) |> 
    mutate(dataset = ds)
})

y_df_nb <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 1/3'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_2_nb.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_pois) |> 
    mutate(dataset = ds)
})

y_df_pois <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 0'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_2_pois.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
fig_cntr <- 1

caption <- str_glue("Fig {fig_cntr}. Scenario 2's simulated incidence reports. Measurement noise from the Poisson (no overdispersion) & Negative Binomial distributions was added to the smooth trajectories obtained from SEIR instances with an exponential-distributed latent period.")
```


```{r, fig.cap = caption}
y_df <- bind_rows(y_df_nb, y_df_pois) |> 
  mutate(highlight = ifelse(dataset == 10, TRUE, FALSE))

plot_y_by_overdispersion(y_df)
```

## Inference

### Poisson data

```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_2_pois"
ll_list    <- vector(mode = "list", length = length(M_j))
```

#### Fitting $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- 1:4
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

    
##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

```{r}
trad_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_2_pois/trad"
create_backup_fldrs(D_i, D_j, 4, M_i, M_j, trad_fldr )
trad_list  <- fit_datasets(D_i, D_j, 4, y_df_pois, info_trad_pois, 
                               trad_fldr, n_chains, samples_per_chain, 
                               n_params = 5, alt_param = FALSE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

In the figure below, we show the results of fitting the **traditional** 
three-unknown parameter model. Here, we can see that the bias caused by the 
misspecification of the infectious period is larger than the one obtained from 
the alternative parameterisation. In other words, it is costlier to misspecify 
the mean generation time than the mean infectious period.

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(trad_list, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0"), 
                            actual_val = c(R0_val, rho_val, 1))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r}
table_coverage(coverage_df, "Table 2. Scenario 2's coverage table. Poisson noise")
```


### Overdispersed data

```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_2_nb"
ll_list    <- vector(mode = "list", length = length(M_j))
```


#### Fitting $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```


```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df, c(0, 5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df, c(0, 5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 3, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```


```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df, c(0, 4))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 3, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df, c(0, 5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 3, x_label = "I[0]", D_i, D_j)
```

#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0", "inv_phi"), 
                            actual_val = c(R0_val, rho_val, 1, inv_phi_ovds))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r}
table_coverage(coverage_df, "Table 3. Scenario 2's coverage table. Overdispersion")
```

# Scenario 3

In comparison to Scenario 2, we increase the reference mean generation time 
($\tau_e$) from **8** to **13**.

## Synthetic data

```{r}
inv_sigma_val <- 8
inv_gamma_val <- 5 # Large
beta_val      <- 0.5 # R0 2.5
ts_val        <- 180
R0_val        <- beta_val * inv_gamma_val
```

```{r}
set.seed(1235)

y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_ovds) |> 
    mutate(dataset = ds)
})

y_df_nb <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 1/3'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_3_nb.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_pois) |> 
    mutate(dataset = ds)
})

y_df_pois <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 0'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_3_pois.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Scenario 3's simulated incidence reports. Measurement noise from the Poisson (no overdispersion) & Negative Binomial distributions was added to the smooth trajectories obtained from SEIR instances with an exponential-distributed latent period.")
```

```{r, fig.cap = caption}
y_df <- bind_rows(y_df_nb, y_df_pois) |> 
  mutate(highlight = ifelse(dataset == 10, TRUE, FALSE))
plot_y_by_overdispersion(y_df)
```


## Inference

### Poisson data

```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_3_pois"
ll_list    <- vector(mode = "list", length = length(M_j))
```

#### Fitting $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```


```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0"), 
                            actual_val = c(R0_val, rho_val, 1))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r}
table_coverage(coverage_df, "Table 4. Scenario 3's coverage table. Poisson noise")
```

### Overdispersed data

```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_3_nb"
ll_list    <- vector(mode = "list", length = length(M_j))
```

#### Fitting $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 4))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 3, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```


```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 3, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 4))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 3, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 4))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 2, 3, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 3, x_label = "I[0]", D_i, D_j)
```

#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0", "inv_phi"), 
                            actual_val = c(R0_val, rho_val, 1, inv_phi_ovds))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r}
table_coverage(coverage_df, "Table 5. Scenario 3's coverage table. Overdispersion")
```

# Scenario 4

In comparison to Scenario 1, we increase $\Re_0$ from **2.5** to **9**.


## Synthetic data

```{r}
inv_sigma_val <- 2
inv_gamma_val <- 2 # low
beta_val      <- 4.5 # R0 9
ts_val        <- 30
R0_val        <- beta_val * inv_gamma_val
```


```{r}
set.seed(18061986)

y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_ovds) |> 
    mutate(dataset = ds)
})

y_df_nb <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 1/3'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_4_nb.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_pois) |> 
    mutate(dataset = ds)
})

y_df_pois <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 0'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_4_pois.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Scenario 4's simulated incidence reports. Measurement noise from the Poisson (no overdispersion) & Negative Binomial distributions was added to the smooth trajectories obtained from SEIR instances with an exponential-distributed latent period.")
```


```{r, fig.cap = caption}
y_df <- bind_rows(y_df_nb, y_df_pois) |> 
  mutate(highlight = ifelse(dataset == 10, TRUE, FALSE))
plot_y_by_overdispersion(y_df)
```


## Inference

### Poisson data

```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_4_pois"
ll_list    <- vector(mode = "list", length = length(M_j))
```

#### Fitting $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois)
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 7, 15, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 0.25))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 7, 15, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois)
```


##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 7, 15, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20 
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(ds, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 0.5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 7, 15, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 5, x_label = "I[0]", D_i, D_j)
```

#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0"), 
                            actual_val = c(R0_val, rho_val, 1))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r S4_pois_cov}
table_coverage(coverage_df, "Table 6. Scenario 4's coverage table. Poisson noise")
```


### Overdispersed data


```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_4_nb"
ll_list    <- vector(mode = "list", length = length(M_j))
```


#### Fitting $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

#### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 40)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1.5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 5, 20, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 10, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1.5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 0, 15, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 8, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1.5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 5, 15, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 8, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidece fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 1.5))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 5, 15, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 4, x_label = "I[0]", D_i, D_j)
```

#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0", "inv_phi"), 
                            actual_val = c(R0_val, rho_val, 1, inv_phi_val))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r}
table_coverage(coverage_df, "Table 7. Scenario 4's coverage table. Overdispersion")
```

# Scenario 5

In comparison to Scenario 4, we increase $\Re_0$ from **9** to **17**.

## Synthetic data

```{r scenario_5_params}
inv_sigma_val <- 2
inv_gamma_val <- 2 # low
beta_val      <- 8.5 # R0 17
ts_val        <- 25
R0_val        <- beta_val * inv_gamma_val
```


```{r}
set.seed(1121)

y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_ovds) |> 
    mutate(dataset = ds)
})

y_df_nb <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 1/3'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_5_nb.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
y_list <- lapply(1:n_ds, function(ds) {
  
  SEIR_measured_incidence(i_val, j_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, inv_phi_pois) |> 
    mutate(dataset = ds)
})

y_df_pois <- do.call(rbind, y_list) |> 
  mutate(disp = "phi^-1 ~'= 0'")

syn_data_file <- "./data/Synthetic/Scenarios/Scenario_5_pois.csv"
save_synthetic_data(syn_data_file, y_df_nb)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Scenario 5's simulated incidence reports. Measurement noise from the Poisson (no overdispersion) & Negative Binomial distributions was added to the smooth trajectories obtained from SEIR instances with an exponential-distributed latent period.")
```

```{r, fig.cap = caption}
y_df <- bind_rows(y_df_nb, y_df_pois) |> 
  mutate(highlight = ifelse(dataset == 10, TRUE, FALSE))
plot_y_by_overdispersion(y_df)
```

## Inference

### Poisson data

```{r}
root_fldr <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_5_pois"
ll_list   <- vector(mode = "list", length = length(M_j))
```

#### Fitting $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 0.25))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 15, 25, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 4, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.caption = caption}
incidence_df <- extract_incidence(posterior_sample, 40)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 0.20))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 10, 25, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 4, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0, 0.4))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 10, 25, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 4, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_pois, info_files_pois, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_pois, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_pois, c(0.8, 1.1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 10, 25, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 4, x_label = "I[0]", D_i, D_j)
```


#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0"), 
                            actual_val = c(R0_val, rho_val, 1))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r}
table_coverage(coverage_df, "Table 8. Scenario 5's coverage table. Poisson")
```

### Overdispersed data

```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/Scenarios/Scenario_5_nb"
ll_list    <- vector(mode = "list", length = length(M_j))
```


#### $D^{11}$

```{r}
D_i    <- 1
D_j    <- 1
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_nb, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 5, 35, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 15, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{12}$

```{r}
D_i    <- 1
D_j    <- 2
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_nb, c(0, 1))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 10, 35, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 8, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{13}$

```{r}
D_i    <- 1
D_j    <- 3
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 1
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 30)
plot_incidence_prediction(incidence_df, y_df_nb, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_nb, c(0, 2))
```

##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```


```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 10, 35, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 8, x_label = "I[0]", D_i, D_j)
```

#### Fitting $D^{14}$

```{r}
D_i    <- 1
D_j    <- 4
n_data <- 20
create_backup_fldrs(D_i, D_j, n_data, M_i, M_j, root_fldr)
```

```{r}
ll_list[[D_j]] <- fit_datasets(D_i, D_j, n_data, y_df_nb, info_files_nbin, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE,
                               inv_sigma = inv_sigma_val,
                               inv_gamma = inv_gamma_val)
```

```{r}
ds               <- 2
posterior_sample <- ll_list[[D_j]][[ds]]

ds_subset        <- c(1, 6, 10, 14)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

##### Incidence fit

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Posterior predictive checks against latent incidence. Dots denote synthetic data while lines indicate simulations from candidates' ODE structure. We configure these structures using samples from the posterior distribution.")
```

```{r, fig.cap = caption}
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df, D_i, D_j, ds)
```

```{r}
fig_cntr <- fig_cntr + 1

caption  <- str_glue("Fig {fig_cntr}. Fit scores by candidate model and type of data. $x$ denotes latent incidence, whereas $y$ indicates observed incidence. Vertical line denotes the mean.")
```

```{r, fig.cap = caption}
plot_hist_mase(posterior_sample, y_df_nb, c(0, 1.5))
```


##### $\Re_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the basic reproduction number by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_R0_comparison_by_dataset(posterior_subset, R0_val, 5, 35, D_i, D_j)
```

##### $I_0$ marginal distribution

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the initial number of infectious individuals by dataset and model candidate. Error bars denote 95\\% credible intervals. The value in the middle of the bars indicates the distance between the marginal posterior's mean and the true value (vertical line).")
```

```{r, fig.cap = caption}
plot_par_estimates_by_dataset(posterior_subset, "I0", 1, x_min = 0,
                              x_max = 12, x_label = "I[0]", D_i, D_j)
```

#### Summary

##### Coverage

```{r}
actual_values <- data.frame(par        = c("R0",   "par_rho", "I0", "inv_phi"), 
                            actual_val = c(R0_val, rho_val,    1  , inv_phi_val))

map_df(ll_list, calculate_coverage, actual_values, inv_gamma_val) |> 
  pivot_wider(names_from = par, values_from = pct_right) -> coverage_df
```

```{r}
table_coverage(coverage_df, "Table 9. Scenario 5's coverage table. Overdispersion")
```