weekly-report.Rmd

---
title: "COVID-19 US Weekly Forecast Summary"
author: <a href="https://covid19forecasthub.org/doc/" target="_blank">The COVID-19 Forecast Hub Team</a> <br><br> <a href="https://covid19forecasthub.org/" target="_blank">https://covid19forecasthub.org/</a>
date: "report generated `r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float:
      collapsed: false
      smooth_scroll: false
params:
  conn: NA
  quiet: false
---
<!-- code to run rmarkdown::render(input="./vignettes/covidHubUtils-overview.Rmd") -->

<!-- Code for adding logo at the top -->

<!-- <script> -->
<!--   $(document).ready(function() { -->
<!--     $('#TOC').parent().prepend('<div id=\"nav_logo\"><a href=\"https://covid19forecasthub.org/\" target=\"_blank\"><img src=\"https://github.com/reichlab/covid19-forecast-hub-web/raw/master/images/forecast-hub-logo_DARKBLUE.png\"></a></div>'); -->
<!--   }); -->
<!-- </script> -->

<!-- <style> -->
<!-- #nav_logo { -->
<!--   width: 100%; -->
<!--   margin-top: 20px; -->
<!-- } -->

<!-- #TOC { -->
<!--   background: url("https://github.com/reichlab/covid19-forecast-hub-web/raw/master/images/forecast-hub-logo_DARKBLUE-20px-padding.png"); -->
<!--   background-size: contain; -->
<!--   padding-top: 80px !important; -->
<!--   background-repeat: no-repeat; -->
<!-- } -->
<!-- </style> -->
<!-- </style> -->

```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = FALSE,
	message = FALSE,
	warning = FALSE
)
library(lubridate)
library(DT)
library(zoltr) ## devtools::install_github("reichlab/zoltr")
library(scico)
library(tidyverse)
library(htmltools)
library(covidHubUtils)
library (zoo)
library(RColorBrewer)
theme_set(theme_bw())
```

```{r zoltar-setup}
zoltar_connection <- NA

# make a new connection if no connection given
if (!inherits(params$conn, "ZoltarConnection")) {
  ## connect to Zoltar
  zoltar_connection <<- new_connection()
  
  # try to connect to Zoltar 5 times; if all fail, fail the script
  num_tries <- 0
  success <- FALSE
  while(num_tries < 5 && !success) {
    tryCatch(
      # try to authenticate
      {
        zoltar_authenticate(
          zoltar_connection,
          Sys.getenv("Z_USERNAME"),
          Sys.getenv("Z_PASSWORD")
        )
        # <<- superassignment: should only do if preceding scope have
        # such a variable!
        # this statement is reached only if authentication is successful
        success <<- TRUE
      },
      # authentication failed! retry
      error = function(c) {
        message(sprintf("Zoltar connection failed! %d retries remaining...", num_tries))
      },
      # add one to number of retries
      finally = function(c) {
        # <<- superassignment: should only do if preceding scope have
        # such a variable!
        num_tries <<- num_tries + 1
      }
    )
  }
  if (!success) {
    stop("Could not make connection to Zoltar after 5 tries")
  }
# use connection if given
} else {
  zoltar_connection <<- params$conn
}


## construct Zoltar query
project_url <- "https://www.zoltardata.com/api/project/44/"
```


```{r get-date-boundaries}
next_saturday <- as.Date(calc_target_week_end_date(today(), horizon= 0))
# use fixed date
# next_saturday <-  as.Date("2022-09-24")

saturday_4_wk_ahead <- next_saturday + 7*3
saturday_4_wk_ahead_txt <- format(saturday_4_wk_ahead, "%B %d, %Y")
saturday_2_wk_ahead <- next_saturday + 7*1
saturday_2_wk_ahead_txt <- format(saturday_2_wk_ahead, "%B %d, %Y")
saturday_1_wk_ahead_txt <- format(next_saturday, "%B %d, %Y")
last_5_saturdays <- next_saturday - 7*c(5:1)
this_monday <- next_saturday - 5
next_2saturday <- next_saturday +7
wk4_monday <- this_monday + 28
wk2_monday <- this_monday + 14
```


# Background
This report provides a brief summary of the weekly ensemble forecast from the <a href="https://covid19forecasthub.org/" target="_blank">COVID-19 Forecast Hub</a> based on forecasts submitted on `r format(this_monday, "%B %d, %Y")`. In collaboration with the US CDC, our team aggregates COVID-19 forecasts from dozens of teams around the globe. Typically on Wednesday of each week, a summary of the week's forecasts from the COVID-19 Forecast Hub appear on the <a href="https://www.cdc.gov/coronavirus/2019-ncov/covid-data/forecasting-us.html" target="_blank">official CDC COVID-19 forecasting page</a>.


```{r nmodels-this-week}
possible_timezeroes <- seq.Date(this_monday, this_monday-6, by="-1 day")
this_week_timezeroes <- timezeros(zoltar_connection, project_url) %>%
  filter(timezero_date %in% possible_timezeroes) %>%
  pull(timezero_date) %>% sort.default()

models_this_week<-load_forecasts(
        # models=c(),
        dates = this_week_timezeroes,
        types = c("point"),
        targets = c("1 day ahead inc hosp"),
        verbose = FALSE)%>%
  pull(model) %>% sort.default()%>% 
  unique()
nmodels_this_week <- length(models_this_week)


```


Every week, teams submit their forecasts to the COVID-19 Forecast Hub. 
This past week, `r nmodels_this_week` models were submitted. 

Each Tuesday, we combine the most recent forecasts from each team into a single "ensemble" forecast of reported COVID-19 hospitalizations at the state and national level. As of February 20, 2023, we are no longer working with case forecasts, and as of March 6, 2023, we are no longer generating ensemble death forecasts.  Older reports include case and death forecasts.  At the moment, we only generate ensemble forecasts for up to four weeks into the future, as  <a href="https://www.pnas.org/doi/full/10.1073/pnas.2113561119" target="_blank">the available evidence</a> suggests that models are less accurate at longer forecast horizons.

An archive of weekly reports from the COVID-19 Forecast Hub can be found at <a href="https://covid19forecasthub.org/doc/reports/" target="_blank">this page</a>. 


```{r count-models}
## how many models in inc_hosp ensemble?
inc_hosp_models <- read_csv(paste0("../../ensemble-metadata/", this_monday, "-inc_hosp-model-weights.csv"))%>%
  select(-locations) %>%
  apply(MARGIN = 2, FUN=function(x) sum(x))
n_inc_hosp_models <- sum(inc_hosp_models>0)

n_unique_models <- length(unique(c(names(inc_hosp_models)[inc_hosp_models>0])))
  
```

```{r aux-data}
locs <- hub_locations %>%
  rename(Population = population)
all_states <-locs$fips[2:58]
```


```{r download-ensemble-data}
inc_hosp_targets <- paste(1:14, "day ahead inc hosp")

# submit query with covidHubUtils

##hospitailization

dat_hosp<-load_forecasts(
        models=c("COVIDhub-ensemble"),
        dates = this_monday, 
        types = c("point", "quantile"),
        targets = paste(inc_hosp_targets),
       location=c("US",all_states),
       verbose = FALSE,
       source = "zoltar") %>%
     rename(fips=location, timezero=forecast_date, class=type, day_ahead=horizon)%>%
     mutate(target=paste(day_ahead,target_variable,sep=" day ahead "))%>%
     select(model, timezero, fips, target, class, quantile, value, day_ahead) %>%
     # create rate variable and week-ahead
     mutate(day_ahead = as.numeric(substr(target, 0,2)),
     ## recreates the target_end_date from GitHub
     target_end_date = as.Date(timezero + day_ahead))    
```


# COVID-19 Forecasts 

## National level {.tabset .tabset-fade}


This week, our ensemble combined forecasts from `r n_unique_models` different models. 

```{r us-summary}

##hospitalization
us_inc_hosp_range <- dat_hosp %>% 
  filter(fips=="US", target %in% inc_hosp_targets, class=="point") %>% 
  pull(value) %>% range() %>%
  round(-2) %>% format(big.mark = ",")

us_inc_hosp_wk_pi_round <- dat_hosp %>% 
  filter(fips=="US", target == "14 day ahead inc hosp", quantile %in% c(0.025, 0.975)) %>% 
  pull(value) %>% sort() %>%
  round(-2) %>% format(big.mark = ",")

us_inc_hosp_wk_pi <- dat_hosp %>% 
  filter(fips=="US", target == "14 day ahead inc hosp", quantile %in% c(0.025, 0.975)) %>% 
  pull(value) %>% sort() %>% 
  format(big.mark = ",")

```
As of September 28, 2021 the ensemble forecast only reports 14 day ahead forecasts for hospitalizations, due to persistent large inaccuracies observed when forecasting beyond that. 

For `r format(wk2_monday, "%B %d, %Y")` COVID-19 daily hospitalizations in the US will be between `r us_inc_hosp_wk_pi_round[1]` and `r us_inc_hosp_wk_pi_round[2]` (95% prediction interval: `r us_inc_hosp_wk_pi[1]` - `r us_inc_hosp_wk_pi[2]`).



<!-- Throughought most of July, models have in general shown broad agreement about the trajectory of the outbreak over the coming weeks. However, the recent surge in cases has left models with quite different interpretations about what the next few weeks hold in terms of how many reported deaths from COVID-19 we will see.  -->
You can explore the full set of models, including their forecasts for past weeks online at the <a href="https://viz.covid19forecasthub.org/" target="_blank">Forecast Hub interactive visualization</a>.

### Hospitalizations

```{r make-US-inc-hospitalization-plot-daily}
quantiles_to_plot <- c(0.025, 0.1, 0.25, 0.75, 0.9, 0.975)

blues <- RColorBrewer::brewer.pal(n=length(quantiles_to_plot)/2+1, "Blues")
inc_hosp_forecast <- dat_hosp %>%
  filter(target %in% inc_hosp_targets)

inc_hosp_truth <- load_truth(
                       truth_source = "HealthData",
                       target_variable = "inc hosp",
                     locations=c("US",all_states))%>%
                    rename(fips = location) %>%
                    left_join(locs, by=c("fips")) %>%
                    filter(fips %in% unique(inc_hosp_forecast$fips), target_end_date <= last_5_saturdays[5])

inc_hosp_all_points <- inc_hosp_truth %>%
  bind_rows(filter(inc_hosp_forecast, class=="point")) %>%
  bind_rows(filter(inc_hosp_truth, target_end_date==last_5_saturdays[5]) %>% mutate(model="COVIDhub-ensemble")) %>%
  mutate(model = relevel(factor(model), ref="Observed Data (HealthData)"))


    
# inc hosp data for code of uncertainty
dummy_inc_hosp <- tibble(
  quantile = quantiles_to_plot, 
  target_end_date=last_5_saturdays[5]) %>%
  right_join(inc_hosp_all_points %>%
      select(-quantile) %>%
      filter(target_end_date == last_5_saturdays[5]))

inc_hosp_quantiles <- inc_hosp_forecast %>%
  dplyr::filter(class=="quantile") %>%
  bind_rows(dummy_inc_hosp) %>%
  dplyr::filter(quantile %in% quantiles_to_plot) %>%
  dplyr::mutate(endpoint_type = ifelse(quantile < 0.5, 'lower', 'upper'),
    alpha = ifelse(endpoint_type == 'lower',
      format(2*quantile, digits=3, nsmall=3),
      format(2*(1-quantile), digits=3, nsmall=3)),
    `Prediction Interval` = fct_rev(paste0((1-as.numeric(alpha))*100, "%"))
  ) %>%
  dplyr::filter(alpha != "1.000") %>%
  dplyr::select(-quantile) %>%
  tidyr::pivot_wider(names_from='endpoint_type', values_from='value')

# daily
ggplot() +
  geom_ribbon(data = inc_hosp_quantiles %>% dplyr::filter(fips=="US"),
    mapping = aes(x = target_end_date,
      ymin=lower, ymax=upper,
      fill=`Prediction Interval`)) +
  geom_line(data=inc_hosp_all_points %>%
      dplyr::filter(fips == "US"),
    mapping = aes(x = target_end_date, y = value, color = model)) +
  geom_point(data=inc_hosp_all_points %>%
      dplyr::filter(fips == "US", !(model=="COVIDhub-ensemble" & target_end_date <= this_monday)),
    mapping = aes(x = target_end_date, y = value, color = model)) +
  scale_fill_manual(values = blues[1:(length(blues)-1)]) +
  scale_color_manual(values = c("black", tail(blues,1))) +
  scale_x_date(name = NULL, date_breaks="6 month", date_labels = "%b %d %Y", date_minor_breaks = "1 month") +
  ylab("incident hospitalizations") +
  labs(title="Daily COVID-19 hospitalization in the US: observed and forecasted",
    caption="source: HealthData (observed data), COVID-19 Forecast Hub (forecasts)") +
  theme(legend.position = c(.05,.95), legend.justification = c(0,1), legend.key = element_rect(colour = "transparent", fill = "white"), legend.background = element_rect(alpha("white", 0.5)), legend.box="horizontal")
```


# Methods & Acknowledgement

This report was reproducibly and dynamically generated using RMarkdown. The code for the report can be found <a href="https://github.com/reichlab/covid19-forecast-hub/tree/master/code/reports" target="_blank">here</a>.


```{r}
# htmltools::includeScript("statcounter.js")
```