epiestim_vignette.Rmd

---
title: "EpiEstim"
output: html_document
date: "2024-09-25"
---

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

# Source functions
```{r source_functions, echo = T, results = 'hide'}
sapply(list.files("../R", full.names = TRUE), source)
```

# Install/load packages

```{r load_packages}
library(devtools)
library(summrt)
library(tidyverse)
library('EpiEstim')
library(ggplot2)
library(lemon)
library(ggtext)
```

# Load in the data

```{r load_data}
all_data <- readRDS("../all_data.RDS")
```

# Prepare the data

Get incidence

```{r get_incidence}
incidence_df <- data.frame(
  dates = all_data$cases$day,
  I = as.vector(all_data$cases$daily_reports)
)
colnames(incidence_df) <- c('dates', 'I')

head(incidence_df)
tail(incidence_df)
dim(incidence_df)
```

Serial interval from data PMF

```{r serial_interval}
si_distr <- as.matrix(all_data$serial$Px)
if (all_data$serial$Day[1] == 1) si_distr <- c(0, si_distr)
si_distr
```

# Run the model

Estimate R daily:

```{r estimate_rt}
getR <- EpiEstim::estimate_R(
  incid = incidence_df,
  method = "non_parametric_si",
  config = EpiEstim::make_config(list(
    si_distr = si_distr,
    t_start = 2:nrow(incidence_df),
    t_end = 2:nrow(incidence_df)
  )),
  backimputation_window = 10
)
```

# Extract the outcomes using `summrt`

```{r extract_outcomes}
output <- summarize_rtestimate(getR)
```

# Plot the data

```{r plot_data_prep}
incubation_shift <- round(weighted.mean(
  x = all_data$incubation$Day,
  w = all_data$incubation$Px
))

reportingdelay_shift <- round(weighted.mean(
  x = all_data$reporting_delay$Day,
  w = all_data$reporting_delay$Px
))
```

Plot data:

```{r plot_data}
# plot_data <- data.frame(
#   package = "EpiEstim",
#   date = all_data$cases$day[getR$R$t_end] - incubation_shift - reportingdelay_shift,
#   Rt_median = getR$R$`Median(R)`,
#   Rt_lb = getR$R$`Quantile.0.025(R)`,
#   Rt_ub = getR$R$`Quantile.0.975(R)`
# )
plot_data <- summarize_rtestimate(getR)
this_plot(plot_data, "EpiEstim")

```