COVID_ensemble_forecasting_anzsc.Rmd

---
title: "Probabilistic ensemble forecasting of Australian COVID-19 cases"
author: "Rob J Hyndman"
date: ISF 2021
fontsize: 14pt
classoption: aspectratio=169
toc: false
output:
  binb::monash:
    fig_height: 4.33
    fig_width: 7
    colortheme: monashwhite
    keep_tex: no
    includes:
      in_header: header.tex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  message = FALSE,
  warning = FALSE,
  cache = TRUE,
  dev.args = list(pointsize = 11)
)
options(digits = 3, width = 88)
library(tidyverse)
library(tsibble)
library(feasts)
library(fable)
library(distributional)
library(gganimate)
source("functions.R")
```

```{r get_local, eval=FALSE, include=FALSE}
localcases <- readr::read_csv("~/git/covid19ensemble/inputs/local_cases_input.csv") %>%
    rename(
      n = count,
      date = date_onset
    ) %>%
    # Adjust count to include cases yet to be detected
    mutate(n = n / detection_probability) %>%
    # Remove last days with prob of detection < 0.5
    filter(detection_probability > 0.5) %>%
    as_tsibble(index = date, key = state) %>%
  select(date, state, n) %>%
  filter(state != "AUS")
saveRDS(localcases,"localcases.rds")
```

```{r read_local, echo=FALSE}
localcases <- readRDS("localcases.rds") 
```

## Australian Health Protection Principal Committee

\begin{block}{}The \textbf{Australian Health Protection Principal Committee} is the key decision-making committee for national health emergencies. It comprises all state and territory Chief Health Officers and is chaired by the Australian Chief Medical Officer.
\end{block}

\begin{alertblock}{COVID-19 forecasting group}
\begin{multicols}{3}\small
\begin{itemize}\tightlist
\item Peter Dawson
\item Nick Golding
\item Rob J Hyndman
\item Dennis Liu
\item James M McCaw
\item Jodie McVernon
\item Pablo \rlap{Montero-Manso}
\item Robert Moss
\item Mitchell \rlap{O'Hara-Wild}
\item David J Price
\item Joshua V Ross
\item Gerry Ryan
\item Freya M Shearer
\item Tobin South
\item Ruarai Tobin
\end{itemize}
\end{multicols}\vspace*{-0.2cm}
\end{alertblock}

## Data sources

* Case-level data of all positive COVID-19 tests: onset and detection times.
* Daily population mobility data from Google, Apple & Facebook
* Weekly non-household contact surveys
* Weekly behavioural surveys
* Daily case numbers from many countries and regions via the Johns Hopkins COVID-19 repository

## Case numbers

```{r, echo=FALSE, fig.height=2.9}
state_colours <- c(
  NSW = "#56b4e9",
  VIC = "#0072b2",
  QLD = "#009e73",
  SA = "#f0e442",
  NT = "#d55e00",
  WA = "#e69f00",
  TAS = "#cc79a7",
  ACT = "#cccccc"
)
localcases %>% 
  autoplot(n) + 
  labs(x="Date of symptom onset") +
  scale_color_manual(values = state_colours)
```

* Recent case numbers are uncertain and incomplete as date of onset is not known until symptoms show and a test is obtained.

## Google mobility data

\vspace*{-0.15cm}\includegraphics[width=15cm,height=7.7cm,keepaspectratio=true,trim=0 0 0 7,clip=true]{google}

\begin{textblock}{7}(8.7,1.4)
\begin{block}{}\small
Percentage change compared to pre-COVID-19 baseline for:
\begin{enumerate}\tightlist
\item[(a)] time at workplace;
\item[(b)] time at retail/recreation;
\item[(c)] time at transit stations.
\end{enumerate}
Vertical lines: physical distancing measures implemented.
\end{block}
\end{textblock}

<!-- ## Facebook mobility data -->

<!-- \full{facebook} -->

<!-- \begin{textblock}{5.3}(10.3,.4) -->
<!-- \begin{block}{}\small -->
<!-- Proportion of Facebook users who “stayed put”, 29 Feb -- 2~Jul~2020. Each line is one LGA. -->
<!-- \end{block} -->
<!-- \end{textblock} -->

## Macrodistancing

\vspace*{-0.15cm}\includegraphics[width=15cm,height=7.7cm,keepaspectratio=true]{macrodistancing_effect}

\begin{textblock}{7}(8.7,1.4)
\begin{block}{}\small
\textbf{Estimated \# non-household contacts per day} based on nationwide weekly surveys (gray) and Google mobility data.\\ Green: public holidays.
\end{block}
\end{textblock}

## Microdistancing

\vspace*{-0.15cm}\includegraphics[width=15cm,height=7.7cm,keepaspectratio=true]{microdistancing_effect}

\begin{textblock}{7}(8.7,1.4)
\begin{block}{}\small
\textbf{Estimated \% keeping 1.5m distance from non-household contacts} based on nationwide weekly surveys (gray).
\end{block}
\end{textblock}

## Global daily cases by region from Johns Hopkins

\begin{alertblock}{}https://github.com/CSSEGISandData/COVID-19
\end{alertblock}

\full{jhu_dashboard}

## Model 1: SEEIIR (Uni Melbourne/Doherty Institute)

* Stochastic susceptible-exposed-infectious-recovered compartmental model that incorporates changes in local transmission potential via a time-varying effective reproduction number.
* Uses mobility and survey data and case onset and detection times.
* Daily counts $\sim$ Negative Binomial. 
* Time in class $\sim$ Gamma. 
* Forecasts obtained using a bootstrap particle filter.

## Model 2: Generative model (Uni Adelaide)

* Uses mobility and survey data and case onset and detection times.
* Three types of infectious individuals: imported, asymptomatic, symptomatic
* Class counts $\sim$ Negative Binomial. 
* Incubation times $\sim$ Gamma. 
* Estimation via Hamilton Monte Carlo
* Forecasts obtained via simulation

## Model 3: Global AR model (Monash)
\fontsize{14}{16}\sf

* Uses Johns Hopkins data from countries and regions with sufficient data.
* Series with obvious anomalies (negative cases and large step changes) removed.
* $n_{t,i}=$ daily cases on day $t$ in country/region $i$ (scaled so all data have same mean and variance).
* $y_{t,i} = \phi_1 y_{t-1,i} + \cdots +\phi_p y_{t-p,i} + \varepsilon_{t,i}$\newline where $y_{t,i} = \log(n_{t,i}+0.5)$ and $\varepsilon_{t,i}\sim N(0,\sigma_i^2)$.
* No stationarity constraints. Common coefficients.
* Current model has $p=24$ (selected to minimize the 7-day-ahead MAE on recent Australian data).

## Forecasting ensemble
\fontsize{14}{17}\sf

* Forecasts obtained from a mixture distribution of the component models.
$$\tilde{p}(y_{t+h}|I_t) = \sum_{k=1}^3 w_{t+h|t,k} p_k(y_{t+h}|I_t)$$
where $p_k(y_{t+h}|I_t)$ is the forecast distribution from model $k$, $I_t$ denotes the data available at time $t$ and the weights $w_{t+h|t,k}>0$ sum to one.
* Also known as "linear pooling"
* Works best when individual models are over-confident and use different data sources.
* We have used equal weights $w_{t+h|t,k}=1/3$.

<!--
## Ensemble forecasts: Victoria

\only<1>{\full{vic_forecasts1}}
\only<2>{\full{vic_forecasts2}}
\only<3>{\full{vic_forecasts3}} -->

## Ensemble forecasts: Victoria

```{r combined_forecasts, eval=FALSE}
# Read weekly samples files from mediaflux and save as rds file
fs::dir_ls("~/mediaflux", glob = "*.csv") %>%
  stringr::str_subset("combined_samples_202") %>%
  purrr::map_dfr(read_csv) %>%
  nest(sample = sim1:sim2000) %>%
  group_by(date, state, .model, forecast_origin) %>%
  mutate(sample = list(unname(unlist(sample)))) %>%
  ungroup() %>%
  saveRDS(file = "samples.rds")
```

```{r read_samples}
samples <- readRDS("samples.rds")
ensemble <- make_ensemble(samples)
```

```{r some_plots, include=FALSE}
vic_ensemble <- ensemble %>% filter(state == "VIC")
origins <- sort(unique(vic_ensemble$forecast_origin))
for (i in seq_along(origins)) {
  p <- vic_ensemble %>%
    filter(forecast_origin == origins[i], date <= origins[i] + 7 * 4) %>%
    mutate(dist = dist_sample(sample)) %>%
    select(-sample) %>%
    as_fable(
      index = date, key = forecast_origin,
      response = "n", distribution = dist
    ) %>%
    autoplot(level = c(50, 60, 70, 80, 90), point_forecast = lst(median)) +
    autolayer(
      filter(localcases, state == "VIC", 
             date >= origins[i] - 7 * 12, date <= origins[i] + 7 * 4),
      n
    ) +
    scale_x_date(
      breaks = seq(as.Date("2020-01-01"), by = "1 month", l = 24),
      minor_breaks = NULL,
      labels = paste(
        rep(c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 2),
        rep(2020:2021, c(12,12))
      )
    ) +
    theme(legend.position = "none") +
    xlab("Date of symptom onset") + ylab("Number of cases")
  pdf(paste0(here::here("figs"), "/ensemble", i, ".pdf"),
    width = 20 / 2.54, height = 10 / 2.54
  )
  print(p)
  crop::dev.off.crop()
}
```

\only<1>{\full{ensemble2}}
\only<2>{\full{ensemble8}}
\only<3>{\full{ensemble11}}
\only<4>{\full{ensemble33}}
\only<5>{\full{ensemble48}}
\only<6>{\full{ensemble88}}
\only<7>{\full{ensemble106}}

<!-- \centerline{\animategraphics[controls,buttonsize=0.3cm,width=14cm]{5}{figs/ensemble}{1}{107}} -->

## Evaluating probabilistic forecasts

\begin{textblock}{9.5}(0.2,1.2)
\begin{alertblock}{}\vspace*{-0.7cm}
\begin{align*}
f_{p,t} &= \text{quantile forecast with prob. $p$ at time $t$.}\\
y_{t} &= \text{observation at time $t$}
\end{align*}
\end{alertblock}\vspace*{-0.3cm}
\uncover<2->{\begin{block}{Quantile score}\vspace*{-0.6cm}
$$
  Q_{p,t} = \begin{cases}
  2(1 - p) \big|y_t - f_{p,t}\big|, & \text{if $y_{t} < f_{p,t}$}\\
  2p \big|y_{t} - f_{p,t}\big|, & \text{if $y_{t} \ge f_{p,t}$} \end{cases}
$$
\end{block}}
\end{textblock}
\begin{textblock}{15}(0.2,5.8)
\uncover<4->{
\begin{itemize}\itemsep=0cm\parskip=0cm
\item Low $Q_{p,t}$ is good
\item Multiplier of 2 often omitted, but useful for interpretation
\item $Q_{p,t}$ like absolute error (weighted to account for likely exceedance)
\item Average $Q_{p,t}$ over $p$ = CRPS (Continuous Ranked Probability Score)
\end{itemize}}
\end{textblock}
\begin{textblock}{6}(10,2)
\only<3->{\animategraphics[loop,autoplay]{10}{COVID_ensemble_forecasting_files/figure-beamer/pinball-}{1}{100}}
\end{textblock}

```{r pinball, eval=FALSE, echo=FALSE, fig.show='animate', interval=1/10, message=FALSE, fig.height=3, fig.width=5, cache=FALSE}
# Turn eval=TRUE to recompute these graphs. They are loaded in the above animategraphics call.
prob <- seq(0.05, 0.95, by = 0.05)
df <- expand.grid(
    error = c(-10, 0, 10),
    p = c(prob, rev(head(prob, -1)[-1]))
  ) %>%
  mutate(
    state = rep(seq(length(p) / 3), rep(3, length(p) / 3)),
    qpt = 2 * p * error * (error > 0) - 2 * (1 - p) * error * (error < 0)
  )
labels <- df %>%
  select(p, state) %>%
  distinct() %>%
  mutate(label = paste0("p = ", sprintf("%.2f", p)))
df %>% ggplot(aes(x = error, y = qpt)) +
  geom_line(aes(group = state), colour = "red") +
  labs(
    x = latex2exp::TeX("Error: $y_t - f_{p,t}$"),
    y = latex2exp::TeX("Q_{p,t}")
  ) +
  geom_label(data = labels, aes(x = 0, y = 17, label = label)) +
  transition_states(state)
```
<!--
## CRPS: Continuous Ranked Probability Score
\vspace*{-0.8cm}\fontsize{13}{14}\sf
\begin{align*}
y_{t} &= \text{observation at time $t$}\\
F_{t}(u) &= \text{Pr}(Y_{t} \le u) = \text{forecast distribution}\\
f_{p,t} &= F^{-1}_t(p) = \text{quantile forecast with prob. $p$}\\
  Q_{p,t} &= \begin{cases}
  2(1 - p) \big|y_t - f_{p,t}\big|, & \text{if $y_{t} < f_{p,t}$}\\
  2p \big|y_{t} - f_{p,t}\big|, & \text{if $y_{t} \ge f_{p,t}$} \end{cases} \\
Y_t \text{ and } Y_t^* &\sim \text{iid with distribution $F_t$.}
\end{align*}

###
\begin{align*}
 \text{CRPS}_t
   & = \int_0^1 Q_{p,t}\, dp \\
   & = \int_{-\infty}^\infty \left[F_t(u) - 1_{y_t \le u}\right]^2 du \\
   & = \textstyle\text{E}|Y_t-y_t| - \frac{1}{2}\text{E}|Y_t-Y_t^*|
\end{align*} -->

## CRPS: Continuous Ranked Probability Score

```{r crps}
crps <- bind_rows(
  ensemble %>% sample_crps(localcases) %>% mutate(Model = "Ensemble"),
  samples %>% filter(.model == "gar") %>% sample_crps(localcases) %>% mutate(Model = "Global AR"),
  samples %>% filter(.model == "moss") %>% sample_crps(localcases) %>% mutate(Model = "SEEIIR"),
  samples %>% filter(.model == "uoa") %>% sample_crps(localcases) %>% mutate(Model = "Generative")
)
```

```{r crps_plot, fig.height=4.6, fig.width=10}
crps %>%
  filter(
    h >= 1, h <= 20,
    state %in% c("NSW", "QLD", "SA", "VIC")
  ) %>%
  ggplot(aes(x = h, y = crps, group = Model, col = Model)) +
  geom_line() +
  facet_wrap(. ~ state, scales="free_y") +
  labs(x = "Forecast horizon (days)", y = "CRPS") +
  scale_color_manual(values = c("#D55E00", "#0072B2","#009E73", "#CC79A7"))
```

\begin{textblock}{14.5}(0.5,8)
\begin{block}{}
For weekly forecasts created from 17 September 2020 to 15 June 2021
\end{block}
\end{textblock}

## What have we learned?

* Diverse models in an ensemble are better than one model, especially when they use different information.
* Understand the data, learn from the data custodians.
* Have a well-organized workflow for data processing, modelling and generation of forecasts, including version control and reproducible scripts.
* Communicating probabilistic forecasts is difficult, but consistent visual design is helpful.

## More information
\fontsize{20}{24}\sf

\href{https://robjhyndman.com}{\faicon{home} robjhyndman.com}

\href{https://twitter.com/robjhyndman}{\faicon{twitter} @robjhyndman}

\href{https://github.com/robjhyndman}{\faicon{github}  @robjhyndman}

\href{mailto:rob.hyndman@monash.edu}{\faicon{envelope}  rob.hyndman@monash.edu}