S1.Rmd

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

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

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

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/synthetic_data.R")
source("./R_scripts/tables.R")
source("./R_scripts/utils.R")
source("./R_scripts/utils_model.R")
```

This electronic supplementary material aims to document in a reproducible manner
the generation of the *instruments* employed to carry out the inference study
presented in the main text. By instruments, we refer to the synthetic incidence
reports and the models used to fit them. Both instruments are the product of a 
common data generating process.


# Synthetic incidence

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

# Population size
N_val         <- params_df |> filter(name == "N") |> pull(value)
# 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
n_ds          <- 20

m_val   <- 1
n_val   <- 1:4
ts_val  <- 60
```

## System (latent) component

The first step for generating synthetic data consists of simulating each of 
the $SE^iI^jR$ instances. As mentioned in the main text, we restrict the number
of instances to eight. Two for the latent period, $i =\{1,3\}$, and four, 
$j =\{1,2,3,4\}$, for the infectious period. Moreover, these instances are 
configured to plausible parameter values and initial conditions.

### Parameter values

Irrespective of the distribution of the latent or the infectious period, all
of the instances are configured to the values below. Notice that the choice of 
these values implies that all instances share an identical $\Re_0$ value (2.5).

<br>

```{r, fig.cap= 'Table 1. Constants'}
par_df <- data.frame(Parameter = c('$$\\beta$$', "$$\\sigma$$", "$$\\gamma$$", "$$\\rho$$", "$$N$$"),
                     Value     = c(beta_val, 1/inv_sigma_val, 1/inv_gamma_val, rho_val, N_val))

table_parameters(par_df)
```

### Initial conditions

Regarding the initial conditions of the eight instances, all states set to zero
except for the number of susceptible individuals ($S$) and infectious 
individuals in the first stage ($I^1$).

<br>

```{r}
inits_df <- data.frame(State = c("$$S$$", "$$I^1$$"),
                       init  = c(N_val - 1, 1))

table_inits(inits_df)
```


### Latent incidence

By latent incidence, we mean the *smooth* incidence predicted by an ODE model. 
The plot below shows the latent incidence predicted by the eight instances of 
the $SE^iI^jR$ framework.

<br>

```{r, warning = FALSE}
i_list <- list(1, 3)
j_list <- list(1:4, 1:4)

map2_dfr(i_list, j_list, function(i, j) {
  SEIR_actual_incidence(i, j, N_val, beta_val, inv_sigma_val, inv_gamma_val, 
                        rho_val, ts_val)
}) -> latent_df
```

```{r, fig.cap = 'Fig 1. Incidence obtained from the ODE instances'}
plot_latent_incidence(latent_df)

ggsave("./plots/Fig_01_Latent_states.pdf", height = 5, width = 5)
```

## Observational component

Subsequently, for each latent incidence, we generate $40$ incidence reports 
using the Negative Binomial distribution. Specifically, we produce *two* sets of
20 reports. The first set's noise level is equivalent to that of the Poisson
distribution (no overdispersion).  Namely, the concentration parameter 
($\phi^{-1}$) is set to zero. We refer to these sets as *high-fidelity*. For the
other sets, we add overdispersion ($\phi^{-1} = 1/3$) and identified them as 
*low-fidelity*.

### High-fidelity (Poisson) $D^{1j}$

```{r}
phi_val <- 0 # Poisson
set.seed(1642)

y_list <- lapply(1:n_ds, function(i) {
  
  SEIR_measured_incidence(m_val, n_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val, phi_val) |> 
    mutate(dataset = i)
})
```

```{r}
caption <- 'Fig 2. Simulated incidence reports. Measurement noise from the Poisson distribution was added to the smooth trajectories obtained from SEIR instances with an exponential-distributed latent period.'
```


```{r, fig.cap = caption}
plot_incidence(y_list, "Poisson")
```



```{r}
syn_data_file <- "./data/Synthetic/SEIR/Case_1.csv"

if(!file.exists(syn_data_file)) {
  write_csv(do.call(rbind, y_list), syn_data_file)
}

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

### Low-fidelity (overdispersed) $D^{1j}$

```{r}
phi_val <- 1/3
set.seed(4192)

y_list <- lapply(1:n_ds, function(i) {
  
  SEIR_measured_incidence(m_val, n_val, N_val, beta_val, inv_sigma_val, 
                          inv_gamma_val, rho_val, ts_val,phi_val) |> 
    mutate(dataset = i)
})
```

```{r}
caption <- 'Fig 3. Simulated incidence reports. Measurement noise from the negative binomial distribution with overdispersion was added to the smooth trajectories obtained from SEIR instances with an exponential-distributed latent period.'
```


```{r,  fig.cap = caption}
plot_incidence(y_list, "Negative Binomial")
```

```{r}
syn_data_file <- "./data/Synthetic/SEIR/Case_2.csv"

if(!file.exists(syn_data_file)) {
  write_csv(do.call(rbind, y_list), syn_data_file)
}

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

```{r}
meas_df <- bind_rows(case_1_df, case_2_df) |> 
  filter(dataset == 5)

g <- plot_meas_incidence_by_overdispersion(meas_df)

ggsave("./plots/Fig_02_Measurements.pdf", plot = g, height = 4, width = 5)
```


### High-fidelity $D^{3j}$

```{r}
i_val   <- 3
j_val   <- 1:4
phi_val <- 0 # Poisson

set.seed(1653)

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, phi_val) |> 
    mutate(dataset = ds)
})
```

```{r}
caption <- 'Fig 4. Simulated incidence reports. Measurement noise from the Poisson distribution was added to the smooth trajectories obtained from SEIR instances with an gamma-distributed latent period.'
```


```{r, fig.cap = caption}
plot_incidence(y_list, "Poisson")
```

```{r}
syn_data_file <- "./data/Synthetic/SEIR/Case_3.csv"

if(!file.exists(syn_data_file)) {
  write_csv(do.call(rbind, y_list), syn_data_file)
}
```

### Low-fidelity $D^{3j}$

```{r}
phi_val <- 1/3
i_val   <- 3
j_val   <- 1:4


set.seed(1928)

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, phi_val) |> 
    mutate(dataset = ds)
})
```

```{r}
caption <- 'Fig 5. Simulated incidence reports. Measurement noise from the negative binomial distribution with overdispersion was added to the smooth trajectories obtained from SEIR instances with an gamma-distributed latent period.'
```


```{r, fig.cap = caption}
plot_incidence(y_list, "Negative Binomial")
```

```{r}
syn_data_file <- "./data/Synthetic/SEIR/Case_4.csv"

if(!file.exists(syn_data_file)) {
  write_csv(do.call(rbind, y_list), syn_data_file)
}
```

# Inference files

One can employ mechanistic models as inference tools to fit data in order to
estimate unknown quantities. In this work, we perform the inference process
through the Statistical software *Stan*. This software requires users to 
specify the instructions for the sampling process in Stan's own language. In the
sections below, we show examples of such language for the three
parameterisations described in the main text.

## Three-unknown parameterisation (traditional)

### Example 1

Here is an example of an $M^{11}$ candidate coupled with a Poisson measurement
model.

```{r}
meas_mdl <- "pois"

i_list <- list(1, 3)
j_list <- list(1:4, 1:4)

stan_fldr <- "./Stan_files/Inference/Three_params/"

info_files <- map2(i_list, j_list, \(M_i, M_j) {
  
  info_files <- create_SEIR_files(M_i, M_j, inv_sigma_val, inv_gamma_val, N_val, 
                                "pois", stan_fldr)
}) |> unlist(recursive = FALSE)

saveRDS(info_files, "./Stan_files/Inference/Three_params/pois/meta_info.rds")
```

```{r}
read_file(info_files[[1]]$stan_filepath) |> cat()
```


### Example 2

Here is an example of an $M^{32}$ candidate coupled with a Negative Binomial
measurement model.

```{r}
meas_mdl <- "nbin"

i_list <- list(1, 3)
j_list <- list(1:4, 1:4)

stan_fldr <- "./Stan_files/Inference/Three_params"

info_files <- map2(i_list, j_list, \(M_i, M_j) {
  
  create_SEIR_files(M_i, M_j, inv_sigma_val, inv_gamma_val, N_val, meas_mdl,
                    stan_fldr)
  
}) |> unlist(recursive = FALSE)

saveRDS(info_files, "./Stan_files/Inference/Three_params/nbin/meta_info.rds")
```

```{r}
read_file(info_files[[7]]$stan_filepath) |> cat()
```

## Four-unknown parameterisation

### Example 1

Here is an example of an $M^{14}$ candidate coupled with a Poisson measurement
model.

```{r}
meas_mdl <- "pois"
i_list   <- list(1)
j_list   <- list(1:4)
stan_fldr <- "./Stan_files/Inference/Four_params"

info_files <- map2(i_list, j_list, \(M_i, M_j) {
  
  info_files <- create_SEIR_files(M_i, M_j, inv_sigma_val, inv_gamma_val, N_val, 
                                "pois", stan_fldr, gamma_unk = TRUE)
}) |> unlist(recursive = FALSE)

saveRDS(info_files, "./Stan_files/Inference/Four_params/pois/meta_info.rds")
```

```{r}
read_file(info_files[[4]]$stan_filepath) |> cat()
```

### Example 2

Here is an example of an $M^{13}$ candidate coupled with a Negative Binomial
measurement model.

```{r}
meas_mdl <- "nbin"
i_list <- list(1)
j_list <- list(1:4)
stan_fldr <- "./Stan_files/Inference/Four_params"

info_files <- map2(i_list, j_list, \(M_i, M_j) {
  
  info_files <- create_SEIR_files(M_i, M_j, inv_sigma_val, inv_gamma_val, N_val, 
                                meas_mdl, stan_fldr, gamma_unk = TRUE)
  
}) |> unlist(recursive = FALSE)

saveRDS(info_files, 
        str_glue("./Stan_files/Inference/Four_params/{meas_mdl}/meta_info.rds"))
```

```{r}
read_file(info_files[[3]]$stan_filepath) |> cat()
```

## Three-unknown parameterisation (Alternative)

### Example 1

Here is an example of an $M^{11}$ candidate coupled with a Poisson measurement
model.

```{r}
meas_mdl <- "pois"
stan_fldr <- "./Stan_files/Inference/Three_params_alt"
i_list   <- list(1)
j_list   <- list(1:4)

info_files <- map2(i_list, j_list, \(M_i, M_j) {
  
  info_files <- create_SEIR_files(M_i, M_j, inv_sigma_val, inv_gamma_val, N_val, 
                                meas_mdl, stan_fldr, gamma_unk = FALSE,
                                alt_param = TRUE)
}) |> unlist(recursive = FALSE)

saveRDS(info_files, 
        str_glue("{stan_fldr}/{meas_mdl}/meta_info.rds"))
```

```{r}
read_file(info_files[[1]]$stan_filepath) |> cat()
```


### Example 2

Here is an example of an $M^{12}$ candidate coupled with a Negative Binomial
measurement model.


```{r}
meas_mdl <- "nbin"
i_list <- list(1)
j_list <- list(1:4)
stan_fldr <- "./Stan_files/Inference/Three_params_alt"

info_files <- map2(i_list, j_list, \(M_i, M_j) {
  
  info_files <- create_SEIR_files(M_i, M_j, inv_sigma_val, inv_gamma_val, N_val, 
                                meas_mdl, stan_fldr, gamma_unk = FALSE,
                                alt_param = TRUE)
}) |> unlist(recursive = FALSE)

saveRDS(info_files, 
        str_glue("{stan_fldr}/{meas_mdl}/meta_info.rds"))
```

```{r}
read_file(info_files[[2]]$stan_filepath) |> cat()
```