S7.Rmd

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

In this appendix, we illustrate the procedure to infer the basic reproduction 
number ($\Re_0$) of influenza during the second wave of the **1918 pandemic** in
Cumberland (Maryland). Specifically, we fit four model candidates from the
**alternative parameterisation** to incidence data. The inference process yields
almost identical $\Re_0$ estimates, regardless of the infectious period
distribution. Furthermore, we compare the results of the alternative 
parameterisation to those of the *traditional* one.


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

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

source("./R_scripts/formulas.R")
source("./R_scripts/Inference.R")
source("./R_scripts/plots.R")
source("./R_scripts/plots_Cumberland.R")
source("./R_scripts/posterior_utils.R")
```

# Incidence data

After enduring a wave of influenza infections during the spring of 1918,
the U.S. Public Health Service organised special surveys in several localities
to determine as accurately as possible the proportion of the population 
infected during the second wave of infections in the autumn of 1918. In the 
figure below, we show the report of new cases detected in the city of 
Cumberland (Maryland) over that period.

```{r, fig.cap = 'Fig 1. Cumberland\'s incidence data'}
Maryland_data <- read_csv("./data/Maryland_incidence.csv", 
                          show_col_types = FALSE)


Cumberland_df <- Maryland_data |> select(Date, Cumberland) |> 
  mutate(time = row_number()) |> 
  rename(y = Cumberland)

plot_epi_bars(Cumberland_df)
```

# Inference

We employ four candidates per parameterisation (traditional and alternative). 
On the one hand, the traditional parameterisation refers to the approach of 
fixing the mean of the epidemiological delays (latent and infectious periods) to
values obtained from the literature, irrespective of their particular 
distribution. On the other hand, the proposed alternative parameterisation 
refers to the special emphasis placed on mean generation time of the **SEIR**, 
while flexibilising the mean and distribution of the epidemiological delays.
Namely, the epidemiological delays can take any mean or shape provided that as a
whole conform to the observed mean generation time. Furthermore, following the 
results shown in **S4** and **S5**, for all candidates, we assume an 
exponentially-distributed latent period ($SE^1I^jR$), where $j = \{1,2,3,4\}$.

## Prior distribution

Prior distributions for both parameterisations correspond to those employed in
*S3*.

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

<!-- ## Three unknown parameters -->

```{r}
root_fldr <- "./Saved_objects/Inference/Cumberland/Three_params"
M_i      <- 1
M_j      <- 1:4


info_files <- readRDS("./Stan_files/Inference/Three_params/nbin/meta_info.rds")[1:length(M_j)]

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

```{r}
trad_df <- fit_single_dataset(info_files, Cumberland_df, root_fldr, n_chains,
                               samples_per_chain, n_params = 3, 
                               alt_param = FALSE, N = 5234,
                               xi = 0.3) |> 
  mutate(tau = mean_generation_time(M_j, 2, 2),
         R0  = par_beta * 2,
         parameterisation = "Traditional")
```


```{r}
ppc_file <- file.path(root_fldr, "./ppc.rds")

if(!file.exists(ppc_file)) {
  map_df(1:4, posterior_pred_checks, 1000, FALSE) -> pred_incidence1 
  saveRDS(pred_incidence1 , ppc_file)
} else {
  pred_incidence1 <- readRDS(ppc_file)
}
```

<!-- ## Three unknown parameters (alternative) -->

```{r}
root_fldr <- "./Saved_objects/Inference/Cumberland/Three_params_alt"
M_i      <- 1
M_j      <- 1:4

info_files <- readRDS("./Stan_files/Inference/Three_params_alt/nbin/meta_info.rds")[1:length(M_j)]

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

```{r}
alt_df <- fit_single_dataset(info_files, Cumberland_df, root_fldr, n_chains,
                               samples_per_chain, n_params = 4, 
                               alt_param = TRUE, N = 5234,
                               xi = 0.3, fixed_tau = 2.85) |> 
  mutate(tau = 2.85,
         R0  = 1 / par_inv_R0,
         parameterisation = "Alternative")
```

```{r}
ppc_file <- file.path(root_fldr, "./ppc.rds")

if(!file.exists(ppc_file)) {
  
  map_df(1:4, posterior_pred_checks, 1000, TRUE) -> pred_incidence2 
  saveRDS(pred_incidence2 , ppc_file)
} else pred_incidence2  <- readRDS(ppc_file)

```


## Posterior distribution

### Incidence fit

```{r}
latent_incidence1 <- extract_incidence(trad_df, 40) |> 
  mutate(parameterisation = "Traditional")

latent_incidence2 <- extract_incidence(alt_df, 40)  |> 
  mutate(parameterisation = "Alternative")

latent_incidence <- bind_rows(latent_incidence1, latent_incidence2)

latent_incidence$parameterisation <- factor(latent_incidence$parameterisation,
                                            levels = c("Traditional", 
                                                       "Alternative"))
```

```{r}
#plot_Cmb_latent_incidence(latent_incidence, Cumberland_df)
```

```{r}
pred_incidence1 <- pred_incidence1 |> mutate(parameterisation = "Traditional")
pred_incidence2 <- pred_incidence2 |> mutate(parameterisation = "Alternative")

pred_incidence <- bind_rows(pred_incidence1, pred_incidence2)

pred_incidence$parameterisation <- factor(pred_incidence$parameterisation,
                                            levels = c("Traditional", 
                                                       "Alternative"))
```

```{r}
caption <- "Fig 2. Posterior predictive checks"
```


```{r, fig.cap = caption}
plot_posterior_pred_checks(pred_incidence, Cumberland_df)
```


```{r}
p1 <- fig_7_main(pred_incidence2, Cumberland_df)
```

### Joint posterior distribution


### Marginal distributions

```{r}
both_df <- bind_rows(trad_df, alt_df)
both_df$parameterisation <- factor(both_df$parameterisation,
                                   levels = c("Traditional", "Alternative"))
```

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

```{r}
caption <- "Fig 3. Estimation of the basic reproduction number by model candidate and parameterisation"
```


```{r, fig.cap = caption}
plot_hist_by_M_j(both_df, "R0", x_label = "R[0]", c(1.8, 2.5), n_bins = 50)
```

```{r}
summary_df <- alt_df |> group_by(M_j) |>  
  summarise(mean_R0     = mean(R0),
            lower_bound = quantile(R0, 0.025),
            upper_bound = quantile(R0, 0.975),
            .groups = "drop")

p2 <- fig_7_inset(summary_df)
```

```{r}
plot_Cumberland <- p1 + inset_element(p2, left = 0.29, bottom = 0.58, right = 0.64, top = 0.99)

ggsave("./plots/Fig_07_Cumberland.pdf", plot_Cumberland , height = 4, 
       width = 5,
       device = cairo_pdf)
```



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

```{r}
caption <- "Fig 4. Estimation of the reporting rate by model candidate and parameterisation"
```

```{r, fig.cap = caption}
plot_hist_by_M_j(both_df, "par_rho", x_label = "rho", c(0, 1), 60)
```

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

```{r}
caption <- "Fig 5. Estimation of the initial number of infectious individuals by model candidate and parameterisation"
```

```{r, fig.cap = caption}
plot_hist_by_M_j(both_df, "I0", x_label = "I(0)", c(0, 8), 30)
```

#### Overdispersion parameter ($\phi^{-1}$)

```{r}
caption <- "Fig 6. Estimation of overdispersion by model candidate and parameterisation"
```

```{r, fig.cap = caption}
plot_hist_by_M_j(both_df, "inv_phi", x_label = "phi^-1", c(0, 1))
```