S1_Text.Rmd

---
title: "S1 Appendix"
output: 
  pdf_document:
    number_sections: true
header-includes:
  - \usepackage{booktabs}
---

A critical requirement for parameter inference on State-Space Models is the
robustness of likelihood estimates. Namely, different runs ought to yield 
similar log-likelihood values from a single set of parameters (model 
configuration). Consequently, in this appendix, we formulate and assess the 
convergence of six measurement model candidates that account for the observed 
incidence. Some models account for the relative contact rate, assuming mobility 
data as a proxy observation. We couple those candidate measurements models with 
two process models. The first process model (PM1) describes the relative contact
rate of the SEI3R within-host profile in terms of Geometric Brownian Motion. The
other process model (PM2) employs the Cox-Ingersoll-Ross formulation to describe
the relative contact rate. Specifically, we perform a sensitivity analysis on 
the accuracy of the likelihood estimates under various particle and integration 
step sizes for each chosen model configuration. Overall, the results allow us to
identify possible model misspecifications and to select an appropriate 
integration step and sample size in a context of finite computational resources.


\tableofcontents 

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

library(doParallel)
library(doRNG)
library(dplyr)
library(extraDistr)
library(kableExtra)
library(imputeTS)
library(lubridate)
library(pomp)
library(purrr)
library(readr)
library(readxl)
library(tidyr)
library(stringr)
library(tictoc)

source("./R_scripts/data.R")
source("./R_scripts/helpers.R")
source("./R_scripts/likelihood_funs.R")
source("./R_scripts/POMP_models.R")
source("./R_scripts/plots.R")

denoms <- 2 ** (2:7)
dts    <- 1 / denoms
```


\newpage

# Data

In this work, we base the inference process on two datasets. The first dataset 
corresponds to the number of COVID-19 cases detected during Ireland's first 
wave, from 29 February 2020 to 17 May 2020. We refer to this dataset as $y^1$
and distinguish between two aggregation levels: daily ($y^1_d$) and weekly 
($y^1_w$). The second dataset ($y^2$) is Apple's driving mobility index which 
quantifies the number of requests made to Apple Maps for directions. These 
indexes are normalised to the value on 28 February 2020. Furthermore, given that
data for May 11-12 is missing, we approximate such values using linear
interpolation. In a similar manner to incidence data, we employ two subsets of
datasets. Accordingly, $y^2_d$ denotes daily measurements, and $y^2_w$ denotes
the measurement of the end of the week.

\hfill


```{r}
data_list <- get_data()
```

```{r}
daily_incidence <- get_irish_data() |> slice(1:79)
wkl_incidence   <- get_wkl_incidence(daily_incidence)

g1              <- daily_epi_trend(daily_incidence, "A) Daily detected cases")
g2              <- weekly_epicurve(wkl_incidence, "B) Weekly detected cases")

dates         <- seq(ymd('2020-02-29'), ymd("2020-05-17"), by = '1 day')
drv_data_obj  <- get_driving_data()
imp           <- drv_data_obj$imputed_data
daily_mob_df  <- data.frame(date = dates, y = imp)
g3            <- plot_daily_mobility(daily_mob_df)
g4            <- plot_wkl_mobility(daily_mob_df)

save_plot <- (g1 | g2) / (g3 | g4)  

dvc      <- "pdf"
fig_path <- str_glue("./paper_plots/Fig_02_Data.{dvc}")

ggsave(fig_path, plot = save_plot, height = 5, width = 5)

dvc      <- "eps"
fig_path <- str_glue("./paper_plots/Fig_02_Data.{dvc}")
ggsave(filename = fig_path, plot = save_plot, device = cairo_ps, height = 5, 
       width = 5)

save_plot
```




\newpage

# Geometric Brownian Motion (DGP1)

## Process model (PM1)

\hfill

\begin{equation}
    \frac{dS}{dt} = - S_t \lambda_t
\end{equation}

\begin{equation}
   \frac{dE}{dt} = S_t \lambda_t - \sigma E_t
\end{equation}

\begin{equation}
   \frac{dP}{dt} = \omega \sigma E_t - \eta P_t
\end{equation}

\begin{equation}
   \frac{dI}{dt} =  \eta P_t - \gamma I_t
\end{equation}

\begin{equation}
   \frac{dA}{dt} =  (1-\omega) \sigma E_t - \kappa A_t
\end{equation}

\begin{equation}
   \frac{dR}{dt} =  \kappa A_t + \gamma I_t
\end{equation}

\begin{equation}
   \lambda_t =  \frac{ \beta_t(I_t + P_t + \mu A_t)}{N_t} 
\end{equation}

\begin{equation}
   \beta_t = \zeta Z_t
\end{equation}

\begin{equation}
   \color{red}
   \frac{dZ}{dt} =  \alpha Z_t dW 
\end{equation}

\begin{equation}
   dW \sim Normal(0, \sqrt{dt})
\end{equation}

## Testing measurement model candidates

### Testing points

We use five point estimates (table below) as probes to assess the reliability of
likelihood estimates. These point estimates correspond to _plausible_ values
obtained in a earlier version of this inference analysis. For the final version,
we assume such values as **given**, and conduct the analysis from the start. 
Considering that these tests are not intended to be exhaustive but rather of 
exploratory nature, any arbitrary point is valid as a probe.

\hfill

```{r}
path      <- "./Saved_objects/SEI3R_GBM/top_10.csv"
test_pars <- read_csv(path)

set.seed(123)
phi_guesses   <- round(rhnorm(10, 0.3), 3)
test_pars$phi <- phi_guesses

pars_list <- transpose(test_pars)


kable_df           <- test_pars
old_names          <- colnames(kable_df)
new_names          <- ifelse(old_names == "P_0", 
                             paste0("$",old_names , "$"),
                             paste0("$\\",old_names , "$"))
                               
colnames(kable_df) <- new_names

kable_df <- kable_df |>  mutate(id = row_number(), .before = everything()) |> 
  slice(1:5)
  

knitr::kable(kable_df, "latex", booktabs = TRUE,
             escape = FALSE, digits = 3)
```

### Summary

Specifically, we formulate six measurement model candidates. We describe
these models using three features:

\hfill

* Whether we fit the model to **daily** or **weekly** data.
* Whether we model incidence measurements using the **Poisson** or the 
**Negative binomial** distribution.
* Whether the model incorporates the mobility data as a proxy measurement for 
the relative effective contact rate. We assume these measurements are normally 
distributed.

\hfill

For each probe (point estimate) and model candidate, we run the particle filter
16 times. We repeat the process for various particle and integration step sizes.
For each batch of 16 runs, we estimate the log-likelihood mean, the 
log-likelihood standard error, and the computational time. 

In a nutshell, the results indicate that when one models daily incidence
measurements with the Poisson distribution, the particle filter does not 
converge as the number of particles increases. This finding suggests model
misspecification under such an assumption.

Furthermore, the accuracy of log-likelihood estimates is not compromised by 
relative large integration steps (e.g., $\frac{1}{8}$), which diminishes the
burden on computational resources. On the contrary, and as expected in Monte 
Carlo simulation, the accuracy relies on the number of samples (particles).

\hfill

```{r}
models_df <- data.frame(
  model_id        = 1:6,
  Frequency       = c("Daily", "Daily", "Daily", "Daily", "Weekly", "Weekly"),
  Incidence_model = c("Pois", "Nbin", "Pois", "Nbin", "Pois", "Pois"), 
  Mobility_model  = c(FALSE, FALSE, TRUE, TRUE, FALSE, TRUE),
  Converges       = c("No", "Yes", "No", "Yes", "Yes", "Yes"))

kable_df <- models_df

colnames(kable_df) <- c("Id", "Frequency", "Incidence", "Mobility", "Converges")

knitr::kable(kable_df, "latex", booktabs = TRUE,
             escape = FALSE, digits = 3)
```

### Candidate 1

This formulation assumes that daily incidence measurements are distributed
according to the Poisson distribution. Also, this structure does not 
incorporate mobility data.

#### Equations

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Pois(C_t) 
\end{equation}

\newpage


#### Convergence test

\hfill

\hfill

This test shows that the likelihood standard error **does NOT** tend to zero as 
the number of particles increases. As a consequence, we cannot proceed with the 
next step in the inference process. That is, using the Iterated Filtering
algorithm for parameter inference inasmuch as this plug-and-play method relies
on robust likelihood estimates. Moreover, this result is also a signal of model misspecification.

\hfill

```{r GBM_1}
model_id  <- 1
folder    <- str_glue("./Saved_objects/SEI3R_GBM/mdl_{model_id}")

tested_model <- models_df[model_id, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("GBM_{model_id}")
seeds        <- c(157050396, 429423692, 339568946, 311085801, 553095796)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_GBM(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, GBM_colour)
```

\newpage

### Candidate 2

This formulation assumes that daily incidence measurements are distributed
according to the Negative binomial distribution. Also, this structure does not 
incorporate mobility data.

#### Equations

\hfill

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Nbin(C_t, \phi^{-1}) 
\end{equation}

#### Time comparison

\hfill

As we will see below, these likelihood estimates converge. For this and other 
models that exhibit convergence, we present a plot comparing computational time
elapsed in order to obtain likelihood estimates for testing point # 1. We do so 
for several scenarios. Here, scenarios refer to the number of particles and the
integration step size.

\hfill


```{r}
plot_ll_time(sens_df |> filter(id == 1), GBM_colour)
```


\newpage

#### Convergence test

\hfill

In this plot, we observe that the likelihood standard error approximates zero as
the number of particles increases.

\hfill

```{r GBM_2}
model_id  <- 2
folder    <- str_glue("./Saved_objects/SEI3R_GBM/mdl_{model_id}")

tested_model <- models_df[model_id, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("GBM_{model_id}")
seeds        <- c(406058540, 384577220, 450611898, 158637695, 445355553)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_GBM(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, GBM_colour)
```

\newpage

### Candidate 3

This formulation assumes that daily incidence measurements are distributed
according to the Poisson distribution. Also, this structure **incorporates**
mobility data.

#### Equations

\hfill

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Pois(C_t) 
\end{equation}

\begin{equation}
  y^2_d \sim Normal(Z_t, \tau) 
\end{equation}

\newpage

#### Convergence test

\hfill

Mobility data does not redress convergence issues in the assumption that daily
incidence measurements are distributed according to the Poisson distribution.

\hfill

```{r GBM_3}
model_id  <- 3
folder    <- str_glue("./Saved_objects/SEI3R_GBM/mdl_{model_id}")

tested_model <- models_df[model_id, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("GBM_{model_id}")
seeds        <- c(749210329, 633239544, 506479957, 675552576, 33000524)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_GBM(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, GBM_colour)
```

\newpage

### Candidate 4

This formulation assumes that daily incidence measurements are distributed
according to the Negative Binomial distribution. Also, this structure 
**incorporates** mobility data.

#### Equations

\hfill

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Nbin(C_t, \phi^{-1}) 
\end{equation}

\begin{equation}
  y^2_d \sim Normal(Z_t, \tau) 
\end{equation}

#### Time comparison

\hfill

\hfill


```{r}
plot_ll_time(sens_df |> filter(id == 1), GBM_colour)
```

\newpage

#### Convergence test

\hfill

\hfill

```{r GBM_4}
model_id  <- 4
folder    <- str_glue("./Saved_objects/SEI3R_GBM/mdl_{model_id}")

tested_model <- models_df[model_id, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("GBM_{model_id}")
seeds        <- c(844589231, 935794189, 132397586, 602352313, 252313981)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_GBM(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, GBM_colour)
```

\newpage

### Candidate 5

This formulation assumes that **weekly** incidence measurements are distributed
according to the Poisson distribution. Also, this structure **does not** 
incorporate mobility data.

#### Equations

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 7)
\end{equation}

\begin{equation}
  y^1_w \sim Pois(C_t) 
\end{equation}

#### Time comparison

\hfill

\hfill

```{r}
plot_ll_time(sens_df |> filter(id == 1), GBM_colour)
```

\newpage

#### Convergence test

\hfill

\hfill

```{r GBM_5}
model_id  <- 5
folder    <- str_glue("./Saved_objects/SEI3R_GBM/mdl_{model_id}")

tested_model <- models_df[model_id, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("GBM_{model_id}")
seeds        <- c(737178508, 153530893, 932012977, 147980594, 153746799)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_GBM(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, GBM_colour)
```

\newpage

### Candidate 6

This formulation assumes that **weekly** incidence measurements are distributed
according to the Poisson distribution. Also, this structure **incorporates**
mobility data.

#### Equations

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 7)
\end{equation}

\begin{equation}
  y^1_w \sim Pois(C_t) 
\end{equation}

\begin{equation}
  y^2_w \sim Normal(Z_t, \tau) 
\end{equation}

#### Time comparison

\hfill

\hfill

```{r}
plot_ll_time(sens_df |> filter(id == 1), GBM_colour)
```

\newpage

#### Convergence test

\hfill

\hfill

```{r GBM_6}
model_id  <- 6
folder    <- str_glue("./Saved_objects/SEI3R_GBM/mdl_{model_id}")

tested_model <- models_df[model_id, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("GBM_{model_id}")
seeds        <- c(155346244,717579098, 627611221, 714870127, 576080780)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_GBM(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, GBM_colour)
```


\newpage

# Cox-Ingersoll-Ross (DGP2)

## Process model (PM2)

\hfill

\begin{equation}
    \frac{dS}{dt} = - S_t \lambda_t
\end{equation}

\begin{equation}
   \frac{dE}{dt} = S_t \lambda_t - \sigma E_t
\end{equation}

\begin{equation}
   \frac{dP}{dt} = \omega \sigma E_t - \eta P_t
\end{equation}

\begin{equation}
   \frac{dI}{dt} =  \eta P_t - \gamma I_t
\end{equation}

\begin{equation}
   \frac{dA}{dt} =  (1-\omega) \sigma E_t - \kappa A_t
\end{equation}

\begin{equation}
   \frac{dR}{dt} =  \kappa A_t + \gamma I_t
\end{equation}

\begin{equation}
   \lambda_t =  \frac{ \beta_t(I_t + P_t + \mu A_t)}{N_t} 
\end{equation}

\begin{equation}
   \beta_t = \zeta Z_t
\end{equation}

\begin{equation}
   \color{red}
   \frac{dZ}{dt} =  \nu(\upsilon - Z_t)  + \sqrt{\alpha}Z_tdW
\end{equation}

\begin{equation}
   dW \sim Normal(0, \sqrt{dt})
\end{equation}

## Testing measurement model candidates

### Testing points

We use five point estimates (table below) as probes to assess the reliability of
likelihood estimates. These point estimates correspond to _plausible_ values
obtained in a earlier version of this inference analysis. For the final version,
we assume such values as **given**, and conduct the analysis from the start. 
Considering that these tests are not intended to be exhaustive but rather of 
exploratory nature, any arbitrary point is valid as a probe.

\hfill


```{r}
path      <- "./Saved_objects/SEI3R_CIR/top_10.csv"
test_pars <- read_csv(path)

set.seed(321)
phi_guesses   <- round(rhnorm(10, 0.3), 3)
test_pars$phi <- phi_guesses

pars_list <- transpose(test_pars)


kable_df           <- test_pars
old_names          <- colnames(kable_df)
new_names          <- ifelse(old_names == "P_0", 
                             paste0("$",old_names , "$"),
                             paste0("$\\",old_names , "$"))
                               
colnames(kable_df) <- new_names

kable_df <- kable_df |>  mutate(id = row_number(), .before = everything()) |> 
  slice(1:5)
  

knitr::kable(kable_df, "latex", booktabs = TRUE,
             escape = FALSE, digits = 3)
```

### Summary

Similarly to the procedure followed in the previous section, we couple the CIR
process model (PM2) with those six measurement model candidates formulated for
the GBM process model (PM1) and repeat the analysis. In doing so, we obtain 
identical insights.

\hfill

```{r CIR_summary}
models_df <- data.frame(
  model_id        = 7:12,
  Frequency       = c("Daily", "Daily", "Daily", "Daily", "Weekly", "Weekly"),
  Incidence_model = c("Pois", "Nbin", "Pois", "Nbin", "Pois", "Pois"), 
  Mobility_model  = c(FALSE, FALSE, TRUE, TRUE, FALSE, TRUE),
  Convergence     = c("No", "Yes", "No", "Yes", "Yes", "Yes"))

kable_df <- models_df

colnames(kable_df) <- c("Id", "Frequency", "Incidence", "Mobility", "Converges")

knitr::kable(kable_df, "latex", booktabs = TRUE,
             escape = FALSE, digits = 3)
```

### Candidate 7

This formulation assumes that daily incidence measurements are distributed
according to the Poisson distribution. Also, this structure does not 
incorporate mobility data.

#### Equations

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Pois(C_t) 
\end{equation}

\newpage


#### Convergence test

\hfill

\hfill

```{r CIR_7}
model_id  <- 7
folder    <- str_glue("./Saved_objects/SEI3R_CIR/mdl_{model_id}")

tested_model <- models_df[model_id - 6, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("CIR_{model_id}")
seeds        <- c(583990633, 778140747, 320782418, 534157346, 904115668)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_CIR(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, CIR_colour)
```

\newpage

### Candidate 8

This formulation assumes that daily incidence measurements are distributed
according to the Negative binomial distribution. Also, this structure does not 
incorporate mobility data.

#### Equations

\hfill

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Nbin(C_t, \phi^{-1}) 
\end{equation}

#### Time comparison

\hfill

\hfill

```{r}
plot_ll_time(sens_df |> filter(id == 1), CIR_colour)
```

\newpage

#### Convergence test

\hfill

\hfill

```{r CIR_8}
model_id  <- 8
folder    <- str_glue("./Saved_objects/SEI3R_CIR/mdl_{model_id}")

tested_model <- models_df[model_id - 6, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("CIR_{model_id}")
seeds        <- c(60473346, 454175526, 714473089, 710821362, 72567013)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_CIR(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, CIR_colour)
```

\newpage

### Candidate 9

This formulation assumes that daily incidence measurements are distributed
according to the Poisson distribution. Also, this structure **incorporates**
mobility data.

#### Equations

\hfill

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Pois(C_t) 
\end{equation}

\begin{equation}
  y^2_d \sim Normal(Z_t, \tau) 
\end{equation}

\newpage

#### Convergence test

\hfill

\hfill

```{r CIR_9}
model_id  <- 9
folder    <- str_glue("./Saved_objects/SEI3R_CIR/mdl_{model_id}")

tested_model <- models_df[model_id - 6, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("CIR_{model_id}")
seeds        <- c(531475785, 704941514, 155949369, 461304304, 993383892)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_CIR(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, CIR_colour)
```

\newpage

### Candidate 10

This formulation assumes that daily incidence measurements are distributed
according to the Negative Binomial distribution. Also, this structure 
**incorporates** mobility data.

#### Equations

\hfill

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 1)
\end{equation}

\begin{equation}
  y^1_d \sim Nbin(C_t, \phi^{-1}) 
\end{equation}

\begin{equation}
  y^2_d \sim Normal(Z_t, \tau) 
\end{equation}

#### Time comparison

\hfill

\hfill

```{r}
plot_ll_time(sens_df |> filter(id == 1), CIR_colour)
```

\newpage

#### Convergence test

\hfill

\hfill

```{r CIR_10}
model_id  <- 10
folder    <- str_glue("./Saved_objects/SEI3R_CIR/mdl_{model_id}")

tested_model <- models_df[model_id - 6, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("CIR_{model_id}")
seeds        <- c(679472152, 354518584, 31755704, 712176479, 17006830)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_CIR(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, CIR_colour)
```

\newpage

### Candidate 11

This formulation assumes that **weekly** incidence measurements are distributed
according to the Poisson distribution. Also, this structure **does not** 
incorporate mobility data.

#### Equations

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 7)
\end{equation}

\begin{equation}
  y^1_w \sim Pois(C_t) 
\end{equation}

#### Time comparison

\hfill

\hfill

```{r}
plot_ll_time(sens_df |> filter(id == 1), CIR_colour)
```

\newpage

#### Convergence test

\hfill

\hfill

```{r CIR_11}
model_id  <- 11
folder    <- str_glue("./Saved_objects/SEI3R_CIR/mdl_{model_id}")

tested_model <- models_df[model_id - 6, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("CIR_{model_id}")
seeds        <- c(8222346, 585928557, 695889124, 788473425, 357515903)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_CIR(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, CIR_colour)
```


\newpage

### Candidate 12

This formulation assumes that **weekly** incidence measurements are distributed
according to the Poisson distribution. Also, this structure **incorporates**
mobility data.

#### Equations

\hfill

\begin{equation}
   \frac{dC}{dt} =  \eta P_t - C_t \delta(t \, mod \, 7)
\end{equation}

\begin{equation}
  y^1_w \sim Pois(C_t) 
\end{equation}

\begin{equation}
  y^2_w \sim Normal(Z_t, \tau) 
\end{equation}

#### Time comparison

\hfill

\hfill

```{r}
plot_ll_time(sens_df |> filter(id == 1), CIR_colour)
```

\newpage

#### Convergence test

\hfill

\hfill

```{r CIR_12}
model_id  <- 12
folder    <- str_glue("./Saved_objects/SEI3R_CIR/mdl_{model_id}")

tested_model <- models_df[model_id - 6, ]

freq         <- tested_model |> pull(Frequency)
inc_mdl      <- tested_model |> pull(Incidence_model)
mob_mdl      <- tested_model |> pull(Mobility_model)
obs_df       <- data_list[[freq]]
mdl_filename <- str_glue("CIR_{model_id}")
seeds        <- c(524909287, 561295907, 741662553, 348746069, 670886625)

map_df(1:5, function(i) {
  
  pars_set <- pars_list[[i]]
  
  map_dfr(dts, function(dt) {
    
    pomp_obj    <- pomp_CIR(inc_mdl, mob_mdl, obs_df, mdl_filename, dt)
    par_obj     <- pomp_obj$pars
    pomp_mdl    <- pomp_obj$mdl
    n_particles <- c(5e3, 1e4, 2e4, 5e4, 1e5)
    fn          <- file.path(folder, str_glue("pf_sensitivity_{i}_dt_{dt}.rds"))
    pars        <- c(unlist(pars_set), par_obj$fixed)
    pomp_mdl    <- pomp_mdl |> pomp(params = pars)
    
    pf_sensitivity(n_particles = n_particles, n_cores = detectCores(), 
                   seed = seeds[[i]], pomp_mdl = pomp_mdl, fn = fn,
                   n_iter = detectCores() * 2) |> 
      mutate(inv_dt = 1/dt)
   }) |> mutate(id = i)
}) -> sens_df
```

```{r, fig.height = 7}
plot_ll_se(sens_df, CIR_colour)
```

\newpage

# Original Computing Environment

```{r}
sessionInfo()
```