S2.Rmd

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

This appendix illustrates the process of fitting various parameterisations
of the $SE^1I^jR$ model ($M^{1j}$) to **high-fidelity $D^{1j}$** incidence 
reports. We mean by *parameterisation* the decision of categorising model
parameters as either *unknown* or *assumed*. For the unknown parameters, we 
construct prior distributions, which will be eventually updated in light of the
data via HMC sampling, resulting in a posterior distribution. On the other hand,
assumed parameters are fixed at their true values.

```{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(rstan)
library(stringr)
library(tictoc)
library(tidyr)

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

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

# Latent period
inv_sigma_val <- params_df |> filter(name == "inv_sigma") |> pull(value)
# Infectious period
inv_gamma_val <- params_df |> filter(name == "inv_gamma") |> pull(value)
# Effective contact rate
beta_val      <- params_df |> filter(name == "beta") |> pull(value)
# Reporting rate
rho_val       <- params_df |> filter(name == "rho") |> pull(value)

R0_val        <- beta_val * inv_gamma_val

y_df          <- read_csv("./data/Synthetic/SEIR/Case_1.csv", 
                          show_col_types = FALSE)
```

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

M_i    <- 1
M_j    <- 1:4
n_data <- 20

R0_bounds <- c(1,4)
```

# Three-unknown parameterisation (traditional)

For candidates from this parameterisation, we assume that two parameters and one
initial condition are unknown: $\beta$, $\rho$, and  $I^1_0$. Each candidate is 
coupled with the appropriate likelihood function. That is, the *Poisson*
distribution.

## Prior distributions {#link1}

The prior distributions shown below apply for all candidates from this 
parameterisation.

### Effective contact rate ($\beta$)

```{r}
fig_cntr <- 1

caption <- str_glue("Fig {fig_cntr}. Histogram from samples obtained from the contact rate's prior distribution")
```


```{r, fig.cap = caption, fig.align='left', fig.height = 3, fig.width = 5}
prior_df <- data.frame(x = rlnorm(1e5))

plot_prior(prior_df, "beta", "beta~'~ lognormal(0,1)'", c(0, 10), 25)
```

### Reporting rate ($\rho$)

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

caption <- str_glue("Fig {fig_cntr}. Histogram from samples obtained from the reporting rate's prior distribution")
```

```{r, fig.cap = caption, fig.align='left', fig.height = 3, fig.width = 5}
prior_df <- data.frame(x = rbeta(1e5, 2, 2))

plot_prior(prior_df, "rho", "rho~'~ beta(2, 2)'", c(0, 1), 15)
```

### Initial number of infectious individuals ($I_0$)

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

caption <- str_glue("Fig {fig_cntr}. Histogram from samples obtained from the prior distribution of the initial number of infectious individuals.")
```

```{r, fig.cap = caption, fig.align='left', fig.height = 3, fig.width = 5}
prior_df <- data.frame(x = rlnorm(1e5))

plot_prior(prior_df, "I[0]", "I[0]~'~ lognormal(0,1)'", c(0, 10), 25)
```

## Posterior distributions

We stratify inference results by the infectious period distribution that
produced the observed incidences. For instance, $D_{11}$ implies that candidate
models fitted incidence reports that stem from the $SE^1I^1R$.

```{r}
M_i        <- 1
n_params   <- 5
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/SEIR/Case_1/Three_params"
info_files <- readRDS("./Stan_files/Inference/Three_params/Pois/meta_info.rds")[1:max(M_j)]

ll_list <- vector(mode = "list", length = length(M_j))
```

### Fitting $D^{11}$

```{r}
D_i <- 1
D_j <- 1
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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params)
```

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

ds_subset        <- c(1, 5, 10, 20)
posterior_subset <- ll_list[[D_j]][ds_subset]
```

#### Incidence fit

The figure below compares actual (points) and simulated (lines) latent incidence
by model candidate.

```{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}
set.seed(123)
incidence_df <- extract_incidence(posterior_sample, 20)
plot_incidence_prediction(incidence_df, y_df, D_i, D_j, ds)
```

To persuade the reader that all models fit the data equally well, we draw on the
mean absolute scaled error (MASE). This quantity is a measure of the accuracy of
forecasts, well-suited for time-series. Therefore, we employ the MASE to compare
each of the 4000 simulated latent incidences against its true counterpart ($x$)
and observed incidence ($y$). By simulated latent incidences, we refer to the 
process of plugging the samples of a posterior distribution into an ODE model to
obtain incidence trajectories via simulation. We summarise the results via
histograms in the plot below. The left column (of panels) contains the 
comparison between simulated and actual latent incidences. The right column of 
panels displays the comparison between simulated latent incidence and the 
observed incidence. Overall, there is no noticeable variation in the histograms
as the fitting model changes (increasing $j$).

```{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)
```

#### Joint posterior distribution

In this plot, we show an example of a joint posterior distribution. 
Specifically, this posterior distribution was derived from fitting the $M^{11}$
candidate to one $D^{11}$ incidence report.

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == D_j) |> 
  select(par_beta, par_rho, I0) 

pairs_posterior(joint_posterior_df, strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions

In this section, we present a comparison between estimated marginal posterior 
distributions (error bars) and actual values (vertical lines). Such 
distributions are the result of fitting  all four candidates to four $D^{11}$ 
incidence reports. The value in the middle of the error bars indicates the
relative error between the actual value and the mean of the marginal posterior
distribution.

##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, 
                              x_label = "I[0]", D_i, D_j)
```

##### Reporting rate ($\rho$) 

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho", 
                              D_i, D_j)
```

### Fitting $D^{12}$

```{r}
D_i <- 1
D_j <- 2
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, info_files, 
                               root_fldr, n_chains, samples_per_chain,
                               n_params = n_params)
```

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

ds_subset        <- c(1, 5, 10, 20)
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}
set.seed(345)
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)
```


```{r}
g2 <- plot_incidence_prediction(incidence_df, 
                                filter(y_df, time %in% c(1, seq(4, 48, 4))), 
                                D_i, D_j, ds, n_col = 1,
                                y_lims = c(0, 600),
                                line_size  = 0.5,
                                point_size = 0.2, 
                                axis_labels_size = 7,
                                axis_x_size = 6,
                                y_lab = NULL,
                                strip_text_size = 7, 
                                sub_size        = 8,
                                caption = FALSE)

 
```

#### Joint posterior distribution

In this plot, we show an example of a joint posterior distribution. 
Specifically, this posterior distribution was derived from fitting the $M^{12}$
candidate to one $D^{12}$ incidence report.

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == D_j) |> 
  select(par_beta, par_rho, I0) 

pairs_posterior(joint_posterior_df, strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions

##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

```{r}
rg2 <- plot_R0_posterior(posterior_sample, R0_val, D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, 
                              x_label = "I[0]", D_i, D_j)
```

##### Reporting rate ($\rho$) 

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho", 
                              D_i, D_j)
```


### Fitting $D^{13}$

```{r}
D_i <- 1
D_j <- 3
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, info_files, 
                               root_fldr, n_chains, samples_per_chain,
                               n_params = n_params)
```

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

ds_subset        <- c(1, 5, 10, 20)
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)
```

#### Joint posterior distribution

In this plot, we show an example of a joint posterior distribution. 
Specifically, this posterior distribution was derived from fitting the $M^{13}$
candidate to one $D^{13}$ incidence report.

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == D_j) |> 
  select(par_beta, par_rho, I0) 

pairs_posterior(joint_posterior_df, strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions

##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, 
                              x_label = "I[0]", D_i, D_j)
```

##### Reporting rate ($\rho$) 

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho", 
                              D_i, D_j)
```

### Fitting $D^{14}$

```{r}
D_i <- 1
D_j <- 4
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, info_files, 
                               root_fldr, n_chains, samples_per_chain,
                               n_params = n_params)
```

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

ds_subset        <- c(1, 5, 10, 20)
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)
```

#### Joint posterior distribution

In this plot, we show an example of a joint posterior distribution. 
Specifically, this posterior distribution was derived from fitting the $M^{14}$
candidate to one $D^{14}$ incidence report.

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_beta, par_rho, I0) 

pairs_posterior(joint_posterior_df,strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions

##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

```{r}
rg4 <- plot_R0_posterior(posterior_sample, R0_val, D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, 
                              x_label = "I[0]", D_i, D_j)
```

##### Reporting rate ($\rho$) 

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho", 
                              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 1. Coverage table")
```



```{r}
ll_df <- map_df(ll_list, function(dataset_fitlist) {
  
  map_df(dataset_fitlist, \(posterior_df) {
    posterior_df |> 
      select(par_beta, par_rho, log_lik, D_i, D_j, M_i, M_j, dataset) |> 
      mutate(tau    = mean_generation_time(M_j, inv_sigma_val, inv_gamma_val),
             R0     = par_beta * inv_gamma_val,
             id     = paste(M_j, dataset, sep = "_"))
  })
})
```


```{r fig3}
f_A <- fig_3A(ll_df, 10, R0_val)

f_B <- plot_incidence_prediction_by_fitting_model(ll_list, y_df, 10, M_i,
                                                  M_j = c(1, 4),
                                                  D_j = c(1, 4))

g <- guide_area()/f_A/f_B + plot_layout(heights = c(0.25, 2.25, 3.25), 
                                        guides = "collect") + 
  plot_annotation(tag_levels = 'A') & 
  theme(plot.tag = element_text(size = 8))

ggsave("./plots/Fig_03_Pois_inf.pdf", plot = g, height = 5.75, width = 5,
       device = cairo_pdf)
```

### MLE criterion

Each bar (column) in the chart below represents the number of times that 
particular model candidate attained the largest likelihood score for a given 
set of incidence reports. Recall that each set of reports comprises 20
time-series. For instance, the first column (left to right) in the first panel
(top-left) indicates that $M^{11}$ outperforms its competitors 12 out of 20 times
(60 %) in fitting $D^{11}$ incidence reports.

```{r}
MLE_df <- MLE_choice(ll_df) |> 
  mutate(D_i = D_i, .before = D_j)
```

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

caption <- str_glue("Fig {fig_cntr}. Score summary")
```

```{r, fig.cap = caption}
plot_MLE_summary(MLE_df)
```


# Four-unknown parameterisation

In this parameterisation, for each candidate, we assume one additional unknown
parameter: the recovery rate ($\gamma$). Furthermore, all candidates are coupled
with the Poisson distribution.

```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/SEIR/Case_1/Four_params"
info_files <- readRDS("./Stan_files/Inference/Four_params/pois/meta_info.rds")[1:max(M_j)]
ll_list    <- vector(mode = "list", length = length(M_j))
M_i        <- 1
n_params   <- 5
```

## Prior distributions

Prior distributions for $\beta$, $\rho$, and $I^1_0$ are identical to those 
described [here](#link1).

### Recovery rate ($\gamma$)

This prior distribution assumes that the magnitude of the mean infectious 
period ($\gamma^{-1}$) is at least one.

```{r}
fig_cntr <- fig_cntr + 1
caption <- str_glue("Fig {fig_cntr}. Histogram from samples obtained from the recovery rate's prior distribution")
```

```{r, fig.cap = caption, fig.align='left', fig.height = 3, fig.width = 5}
prior_df <- data.frame(x = rbeta(1e5, 2, 2))

plot_prior(prior_df, "gamma", "gamma~'~ beta(2, 2)'", c(0, 1), 15)
```

## Posterior distributions

### 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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params)
```

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

ds_subset        <- c(2, 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)
```


#### Joint posterior distribution {#link2}

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_beta, par_rho, I0, par_gamma) 

pairs_posterior(joint_posterior_df,strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

##### Reporting rate ($\rho$) 

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho", D_i, 
                              D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, x_label = "I[0]", D_i, D_j)
```

##### Recovery rate ($\gamma$)

For this parameter, we opt for a different visualisation. Specifically, for 
each model candidate fitting a particular $D^{11}$ incidence report, we 
compare the marginal posterior distribution against the prior
distribution via histograms. These plots show that the data does little to 
inform the prior distribution as it can be seen that both distributions agree
to a large extent along the parameter space. For this reason, we claim that
the four-unknown parameterisation is **unidentifiable**. This statement is also
supported by the ridge created by $\beta$ and $\gamma$ observed [here](#link2).

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Comparison between the recovery rate's prior distribution (grey histograms) against four marginal posterior distributions obtained from each candidate (coloured histograms).")
```

```{r, fig.cap = caption}
prior_df <- data.frame(sims = rbeta(4e3, 2, 2))
prior_posterior_comparison(posterior_sample, prior_df, "par_gamma")
```

### 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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params)
```

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

ds_subset        <- c(2, 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)
```

#### Joint posterior distribution

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_beta, par_rho, I0, par_gamma) 

pairs_posterior(joint_posterior_df, strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 6
                              , D_i, D_j)
```

##### Reporting rate ($\rho$)

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho",
                              D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, x_label = "I[0]", D_i, D_j)
```

##### Recovery rate ($\gamma$)

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Comparison between the recovery rate's prior distribution (grey histograms) against four marginal posterior distributions obtained from each candidate (coloured histograms).")
```

```{r, fig.cap = caption}
prior_df <- data.frame(sims = rbeta(4e3, 2, 2))
prior_posterior_comparison(posterior_sample, prior_df, "par_gamma")
```

### 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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params)
```

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

ds_subset        <- c(2, 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)
```

#### Joint posterior distribution

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_beta, par_rho, I0, par_gamma) 

pairs_posterior(joint_posterior_df, strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 8, D_i, 
                              D_j)
```

##### Reporting rate ($\rho$)

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho",
                              D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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)
```

##### Recovery rate ($\gamma$)

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Comparison between the recovery rate's prior distribution (grey histograms) against four marginal posterior distributions obtained from each candidate (coloured histograms).")
```

```{r, fig.cap = caption}
prior_df <- data.frame(sims = rbeta(4e3, 2, 2))
prior_posterior_comparison(posterior_sample, prior_df, "par_gamma")
```

### Fitting $D^{14}$

```{r}
D_i <- 1
D_j <- 4
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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params)
```

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

ds_subset        <- c(2, 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}
set.seed(1158)
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)
```


#### Joint posterior distribution

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_beta, par_rho, I0, par_gamma) 

pairs_posterior(joint_posterior_df, strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```


#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 6, D_i,
                              D_j)
```

##### Reporting rate ($\rho$)

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Estimates of the reporting rate 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, "par_rho", rho_val, 
                              x_min = 0.5, x_max = 1, x_label = "rho",
                              D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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)
```

##### Recovery rate ($\gamma$)

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Comparison between the recovery rate's prior distribution (grey histograms) against four marginal posterior distributions obtained from each candidate (coloured histograms).")
```

```{r, fig.cap = caption}
prior_df <- data.frame(sims = rbeta(4e3, 2, 2))
prior_posterior_comparison(posterior_sample, prior_df, "par_gamma")
```

## Summary

### Coverage


```{r}
actual_values <- data.frame(par = c("R0", "par_rho", "I0", "inv_gamma"), 
                            actual_val = c(R0_val, rho_val, 1, inv_gamma_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 2. Coverage table")
```

### $\Re_0$ vs $\tau$

```{r}
ll_df <- map_df(ll_list, function(dataset_fitlist) {
  
  map_df(dataset_fitlist, \(posterior_df) {
    posterior_df |> 
      select(par_beta, par_gamma, par_rho, log_lik, D_i, D_j, M_i, M_j, dataset) |> 
      mutate(tau    = mean_generation_time(M_j, inv_sigma_val, 1 / par_gamma),
             R0     = par_beta / par_gamma,
             id     = paste(M_j, dataset, sep = "_"))
  })
})
```

```{r}
fig_cntr <- fig_cntr + 1
caption  <- str_glue("Fig {fig_cntr}. Scatterplot of the interaction between the mean generation time and the basic reproduction number.")
```

```{r, fig.cap = caption, fig.height = 6}
ll_subset <- ll_df |> filter(dataset %in% c(1, 10, 15))

plot_scatterplot_R0_vs_tau(ll_subset, c(0,15), c(0, 10))
```


```{r}
f_B1 <- fig_4B(ll_df, 20, c(0,10), c(0, 6))
```

```{r}
f_B2 <- fig_4C(ll_df, 20, c(0, 10), c(0, 6))
```


```{r}
design <- "
  111##
  11122
  11122
  11122
  111##"

f_B <- (f_B1|f_B2) + plot_layout(design = design)
```

```{r}
f_A <- fig_4A(ll_df, 20, R0_val) 
```


```{r}
g <- guide_area()/f_A/f_B + plot_layout(heights = c(0.25, 1.75, 4.5),
                           guides = "collect") + 
  plot_annotation(tag_levels = 'A') & 
  theme(plot.tag     = element_text(size = 8))

ggsave("./plots/Fig_04_Four_unknowns.pdf", plot = g, height = 6.5, width = 5, 
      device = cairo_pdf)
```



# Three-unknown parameterisation (Alternative)

The *alternative parameterisation* refers to the algebraic manipulation of the
$SE^iI^jR$ framework so as to obtain a set equations wherein the basic 
reproduction number ($\Re_0$) and the mean generation time ($\tau$) are 
explicit parameters of the model. In doing so, $\beta$ and $\gamma$ become
functions of other variables rather than parameters. In these parameterisation,
we assume the basic reproduction number's inverse, the reporting
rate ($\rho$) and initial number of infectious individuals ($I^1_0$) as 
unknowns. The remaining parameters are fixed to their true values.

\begin{equation}
  \begin{aligned}
    \beta       &= \frac{2 j (\tau - \sigma^{-1})}{\Re_0^{-1} (j + 1)}\\
    \gamma      &= \beta \Re_0^{-1}
  \end{aligned}
\end{equation}

\begin{equation}
  \begin{aligned}
    \dot{S_t}   &= -\frac{\beta S_t \sum^j_{k=1} I^k_t}{N}\\
    \dot{E^1_t} &= \frac{\beta S_t \sum^j_{k=1} I^k_t}{N} - i\sigma E^1_t \\
    \dot{E^2_t} &= i\sigma E^1_t - i\sigma E^2_t \\
    \vdots \\
    \dot{E^i_t} &= i \sigma E^{i - 1}_t - i \sigma E^i_t \\
    \dot{I^1_t} &= i \sigma E^{i}_t - \gamma I^1_t \\
    \dot{I^2_t} &= j \gamma I^{1}_t - \gamma I^2_t  \\
    \vdots \\
    \dot{I^j_t} &= j \gamma I^{j - 1}_t - j \gamma I^{j}_t \\
    \dot{R_t}   &= j \gamma I^j_t
  \end{aligned}
\end{equation}


```{r}
root_fldr  <- "./Saved_objects/Inference/Synthetic_data/SEIR/Case_1/Three_params_alt"
info_files <- readRDS("./Stan_files/Inference/Three_params_alt/pois/meta_info.rds")[1:max(M_j)]
ll_list    <- vector(mode = "list", length = length(M_j))
M_i        <- 1
n_params   <- 5
```

## Prior distributions {#link1}

The prior distributions shown below apply for all candidates from this 
parameterisation. Prior distributions for $\rho$, and $I^1_0$ are identical to 
those described [here](#link1).


### The inverse of the basic reproduction number ($\Re_0^{-1}$)

Since model candidates are fitting outbreak-like incidence data, it is warranted
to assume that the basic reproduction number is higher than one ($\Re_0 > 1$).
This observation implies that the magnitude of its inverse ($\Re_0^{-1}$) 
between 0 and 1. Consequently, we formulate a prior that is consistent with
such a constraint.

```{r}
fig_cntr <- fig_cntr + 1
caption <- str_glue("Fig {fig_cntr}. Histogram from samples obtained from the prior distribution of the basic reproduction number's inverse.")
```

```{r, fig.cap = caption, fig.align ='left', fig.height = 3, fig.width = 5}
prior_df <- data.frame(x = rbeta(1e5, 2, 2))

plot_prior(prior_df, "\u211c^-1", "\u211c^-1~'~ beta(2, 2)'", c(0, 1), 15)
```

## Posterior distributions

### 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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE)
```

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

ds_subset        <- c(3, 6, 12, 18)
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)
```

#### Joint posterior distribution

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_inv_R0, par_rho, I0) 

pairs_posterior(joint_posterior_df,strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, 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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE)
```

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

ds_subset        <- c(3, 6, 12, 18)
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)
```

#### Joint posterior distribution

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_inv_R0, par_rho, I0) 

pairs_posterior(joint_posterior_df,strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, 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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE)
```

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

ds_subset        <- c(3, 6, 12, 18)
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)
```

#### Joint posterior distribution

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_inv_R0, par_rho, I0) 

pairs_posterior(joint_posterior_df,strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, 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, info_files, 
                               root_fldr, n_chains, samples_per_chain, 
                               n_params = n_params, alt_param = TRUE)
```

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

ds_subset        <- c(3, 6, 12, 18)
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)
```

#### Joint posterior distribution

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

caption  <- str_glue("Fig {fig_cntr}. Posterior distribution. The diagonal shows the posterior marginal distributions (histograms). In the lower triangular part, the joint posterior distribution of each possible\ncombination of two parameters is displayed (heatmaps). The upper triangular part shows the correlation among parameters.")
```

```{r, fig.cap = caption}
joint_posterior_df <- posterior_sample |> filter(M_j == !!D_j) |> 
  select(par_inv_R0, par_rho, I0) 

pairs_posterior(joint_posterior_df,strip_text = 10, axis_text_size = 7,
                M_j = D_j)
```

#### Marginal posterior distributions


##### Basic reproduction number ($\Re_0$)

```{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, R0_bounds[[1]], 
                              R0_bounds[[2]], D_i, D_j)
```

##### Initial number of infectious individuals ($I_0$)

```{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 = 2, x_label = "I[0]", D_i, D_j)
```

## Summary

```{r}
ll_df <- map_df(ll_list, function(dataset_fitlist) {
  
  map_df(dataset_fitlist, \(posterior_df) {
    posterior_df |> 
      select(par_inv_R0, I0, par_rho, log_lik, D_j, M_i, M_j, dataset) |> 
      mutate(tau    = mean_generation_time(D_j, inv_sigma_val, inv_gamma_val),
             R0     = 1 / par_inv_R0,
             id     = paste(M_j, dataset, sep = "_"))
  })
})
```


```{r}
g <- fig_5(ll_df, 20, R0_val) 
```

```{r}
ggsave("./plots/Fig_05_Three_unknowns_alt.pdf", plot = g, height = 2.75, 
       width = 5, device = cairo_pdf)
```



### 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 3. Coverage table")
```

<br>