ViroReportR

The goal of ViroReportR is to provide a toolbox to conveniently generate short-term forecasts (with accompanied diagnostics) for viral respiratory diseases.

Installation

ViroReportR depends on the latest version of the EpiEstim package (2.4). Thus, this version of the package must be installed from GitHub prior to installing the ViroReportR package using:

# install.packages("devtools")
install.packages('EpiEstim', repos = c('https://mrc-ide.r-universe.dev', 'https://cloud.r-project.org'))

You can then install the development version of ViroReportR from GitHub with:

devtools::install_github("BCCDC-PHSA/ViroReportR")

Quick Start

ViroReportR can be used to generate short-term forecasts with accompanied diagnostics in a few lines of code. We go through an example here where the EpiEstim backend is used to generate forecasts of Influenza-A, RSV and SARS-CoV-2.

library(ViroReportR)
library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr     1.1.4     ✔ readr     2.1.5
#> ✔ forcats   1.0.0     ✔ stringr   1.6.0
#> ✔ ggplot2   4.0.0     ✔ tibble    3.2.1
#> ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
#> ✔ purrr     1.0.2     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(here)
#> here() starts at C:/Users/rebeca.falcao/ViroReportR
library(DT)
library(purrr)
library(kableExtra)
#> 
#> Attaching package: 'kableExtra'
#> 
#> The following object is masked from 'package:dplyr':
#> 
#>     group_rows

We will use simulate_data to simulate date for Influenza A, RSV, and SARS-CoV-2, which is included with the ViroReportR package. We then pivot the simulated data to transform into a dataset with three columns: date, disease_type and confirm in accordance to format accepted by the model fitting functions.

diseases<-c("flu_a", "rsv", "sars_cov2")
data <- simulate_data(days=365, #days spanning simulation
                      peaks = c("flu_a"=90,"rsv"=110,"sars_cov2"=160), #peak day for each disease
                      amplitudes=c("flu_a"=50,"rsv"=40,"sars_cov2"=20), #amplitude of peak for each disease
                      scales = c("flu_a"=-0.004,"rsv"=-0.005,"sars_cov2"=-0.001), # spread of peak for each disease
                      time_offset = 0, #number of days to offset start of simulation
                      noise_sd = 5, #noise term
                      start_date = "2024-01-07" #starting day (Sunday)
                      )

data$date <- lubridate::ymd(data$date)

vri_data_list <- map2(rep(list(data), length(diseases)),
  diseases,~ get_aggregated_data(.x, "date", .y)) %>% set_names(diseases)

head(vri_data_list)
#> $flu_a
#> # A tibble: 366 × 2
#>    date       confirm
#>    <date>       <dbl>
#>  1 2024-01-07       0
#>  2 2024-01-08       8
#>  3 2024-01-09       0
#>  4 2024-01-10       5
#>  5 2024-01-11       0
#>  6 2024-01-12       8
#>  7 2024-01-13       9
#>  8 2024-01-14       0
#>  9 2024-01-15       4
#> 10 2024-01-16       0
#> # ℹ 356 more rows
#> 
#> $rsv
#> # A tibble: 366 × 2
#>    date       confirm
#>    <date>       <dbl>
#>  1 2024-01-07       0
#>  2 2024-01-08       0
#>  3 2024-01-09       0
#>  4 2024-01-10       0
#>  5 2024-01-11       5
#>  6 2024-01-12       0
#>  7 2024-01-13       0
#>  8 2024-01-14       1
#>  9 2024-01-15       0
#> 10 2024-01-16       0
#> # ℹ 356 more rows
#> 
#> $sars_cov2
#> # A tibble: 366 × 2
#>    date       confirm
#>    <date>       <dbl>
#>  1 2024-01-07       9
#>  2 2024-01-08       0
#>  3 2024-01-09       0
#>  4 2024-01-10       0
#>  5 2024-01-11       0
#>  6 2024-01-12       0
#>  7 2024-01-13       0
#>  8 2024-01-14       3
#>  9 2024-01-15       0
#> 10 2024-01-16       0
#> # ℹ 356 more rows

Model fitting and forecasting

The code below estimates the reproduction number using EpiEstim through the generate_forecast() function for each disease type via the purrr::map2() function. The generate_forecast() function prepares the data, estimates the reproduction number, and produces short-term forecasts of daily confirmed cases for an n_days forecast horizon. The other current choice for the forecasting algorithm is EpiFilter (WIP).

# parameters set-up 
start_date <- min(data$date) + 13
n_days <- 14 # number of days ahead to forecast (n_days)
smooth <- FALSE # logical indicating whether smoothing should be applied in the forecast

forecasts_results <- tibble(
  vri_data_list,
  forecasts = map2(
    vri_data_list,
    diseases,
    ~ generate_forecast(
      data = .x,
      smooth_data = smooth,
      type = .y,
      n_days = n_days,
      start_date = start_date
    )
  )
)


names(forecasts_results$forecasts) <- diseases
names(forecasts_results$vri_data_list) <- diseases

Plotting results

The code below plots the forecasts results and the estimated $R_t$ for each disease. To plot $R_t$, the code below uses plot_rt function included in the package.

for (vri in diseases) {
  forecast_plot <- ggplot() +
  geom_ribbon(
    data = forecasts_results$forecasts[[vri]][["forecast_res_quantiles"]],
    aes(x = date, ymin = p10, ymax = p90, fill = "Prediction Interval"),
    alpha = 0.4
  ) +
  geom_line(
    data = forecasts_results$forecasts[[vri]][["forecast_res_quantiles"]],
    aes(x = date, y = p50, color = "Median Forecast"),
    linewidth = 1
  ) +
  geom_point(
    data = forecasts_results$vri_data_list[[vri]],
    aes(x = date, y = confirm, shape = "Positive Tests"),
    size = 2,
    color = "black"
  )  +
  
  labs(
    title = paste0(n_days, "-Day Forecast of Daily ", vri, " Positive Tests"),
    x = "Date", y = "Predicted Tests",
    fill = "", color = "", shape = ""
  ) +
  scale_fill_manual(values = c("Prediction Interval" = "skyblue")) +
  scale_color_manual(values = c("Median Forecast" = "blue")) +
  scale_shape_manual(values = c("Positive Tests" = 16)) +
  theme_minimal() +
  theme(legend.position = "top", legend.direction = "horizontal")
  
  # create Rt plot
  rt_plot <- ViroReportR:::plot_rt(forecasts_results$forecasts[[vri]])
  print(forecast_plot)
  print(rt_plot)
}

Forecast report

Finally, the ViroReportR package can generate an automated report for the current season across all supported respiratory viruses (Influenza A, RSV, and SARS-CoV-2) using the generate_forecast_report() function. This function renders an HTML report summarizing model outputs and forecasts.

To use it, provide an input file (input_file) containing the required data with three columns—date, disease_type, and confirm—and specify an output directory (output_directory) where the report will be saved.

# rendering forecast report
df <- imap_dfr(vri_data_list, ~ .x %>% mutate(disease_type = .y))
input_file<-"simulated_data.csv"
write.csv(df, input_file, row.names = FALSE)
generate_forecast_report(
  input_data_dir =  input_file, # input filepath
  output_dir = output_directory, # output directory
  n_days = 14, # number of days to forecast
  validate_window_size = 7, # number of days between each validation window
  smooth = FALSE, # logical indicating whether smoothing should be applied in the forecast
)