cityforecasts

R scripts to generate city-level forecasts

Summary

This repository contains R code to generate city-level forecasts for submission to the flu-metrocast hub. All models are exploratory/preliminary, though we will regularly update this document to describe the latest mathematical model used in the submission.

All outputs submitted to the Hub will be archived in this repository, along with additional model metadata (such as the model definition associated with a submission and details on any additional data sources used or decisions made in the submission process). If significant changes to the model are made during submission, we will rename the model in the submission file.

Since the Hub now solicits both local level and aggregate (typically state-level) forecasts of the specified target (typically the percent of ED visits due to flu), we will:

• run each set of local jurisdictions within an aggregate location jointly
• run each aggregate location (typically state) independently using the aggregate data

A table of all the locations being solicited is available at the flu-metrocast hub GitHub.

Since all forecast targets are percents, we will use the same latent and observation model for all locations (however, if we have count data this can easily be incorporated by using a Poisson observation model).

Workflow

Because all data is available publicly, the forecasts generated should be completely reproducible for any given forecast date. We use the mvgam package, which is a an R package that leverages both mgcv and brms formula interface to fit Bayesian Dynamic Generalized Additive Models (GAMs). These packages use metaprogramming to produce Stan files, and we also include the Stan code generated by the package.

To produce forecasts each week we follow the following workflow:

1. Modify the forecast_date in targets/create_config_targets.R. If using for a real-time submission, set real_time <- TRUE, filepath_forecasts <- cityforecasts.
2. Run targets::tar_make(). This should take about ~10 minutes to run on all 70 locations (some of which are estimated jointly). It should populate output/{filepath_forecasts}/{forecast_date}.
3. Open a PR into the flu-metrocast hub GitHub GitHub under our team name epiforecasts.

Steps 1-2 are automated with GitHub Actions, set on a schedule to run at 3:00 PM EST (4:00 PM EDT) every Wednesday. They push the forecasts to a forecast/{forecast_date} branch, publish a report for review on a public GitHub pages site at https://epiforecasts.io/cityforecasts/{forecast_date}, and open a PR into main linking to the report, with a review requested from @kaitejohnson.

Step 3 is automated too, but only runs once a person has reviewed the forecasts. Merging a Forecasts for {forecast_date} PR into main triggers the Submit to hub workflow. It commits the forecast to a new branch off the hub's main, pushes it to the kaitejohnson/flu-metrocast fork and opens a PR into the hub. To submit (or resubmit) a forecast that is already on main, run the workflow manually from the Actions tab with its forecast date, or run FORECAST_DATE={forecast_date} bin/submit_to_hub.sh locally (DRY_RUN=1 commits without pushing).

Reviewing past forecasts

Reports from every automated run are kept on GitHub Pages, so forecasts from past forecast dates can be reviewed at any time:

• The index page, https://epiforecasts.io/cityforecasts/, lists every report by forecast date, newest first.
• Each report is at https://epiforecasts.io/cityforecasts/{forecast_date}/, where {forecast_date} is the Wednesday forecast date in YYYY-MM-DD format (e.g. https://epiforecasts.io/cityforecasts/2026-10-14/).
• Each weekly Forecasts for {forecast_date} PR also links to its report.

The reports are stored on the gh-pages branch. A rerun for the same forecast date replaces that date's report.

Model definition

The below describes the preliminary model used:

Observation model

For the forecasts of flu hospital admissions (counts), we assume a Poisson observation process

$$ y_{l,t} \sim Poisson(exp(x_{l,t})) $$

For the forecasts of the percent of ED visits due to flu, we assume a lognormal observation process on the proportion of ED visits due to flu:

$$\begin{align} p_{l,t} = y_{l,t} \times 100 \\\ log(y_{l,t}) \sim Normal(x_{l,t}, \sigma_{obs,l}) \end{align}$$

The lognormal assumes constant relative (multiplicative) observation noise, so that changes at low proportions (e.g. early-season upticks) are attributed to the latent state rather than treated as noise, as a Beta observation model with constant precision does. Because the reported proportions are rounded to 0.0001, observed zeros (true value < 0.00005) are replaced with 0.00005 before fitting.

Latent state-space model: Dynamic hierarchical GAM with independent autoregression

We model the latent state with a hierarchical GAM component to capture shared seasonality and a univariate autoregressive component to capture deviations from the seasonal pattern within each location.

$$\begin{align} x_{l,t} \sim Normal(\mu_{l,t} + \delta_{l} (x_{l,t-1} - \mu_{l,t-1}), \sigma_l)\\\ \mu_{l,t} = \beta_l + f_{global,t}(weekofyear) + f_{l,t}(weekofyear) \\\ \beta_l \sim Normal(\beta_{global}, \sigma_{count}) \\\ \sigma_{count} \sim exp(0.33) \\\ \delta_l \sim Normal(0.5, 0.25) T[0,1] \\\ \sigma_l \sim exp(1) \\\ \end{align}$$

When a location is fit on its own (e.g. aggregate/state locations), there is no hierarchy: $\mu_{t} = \beta + f(weekofyear)$ and the prior below is placed directly on $\beta$.

For the forecasts of flu hospital admissions (counts), $\beta_{global}$ represents the intercept on the log count scale, so we centre its prior on the log of the typical count, averaging the log of each location's median count:

$$ \beta_{global} \sim Normal(\frac{1}{L}\sum_{l=1}^L log(max(median_t(y_{l,t}), 1)), 2) $$

For the forecasts of percent of ED visits due to flu, $\beta_{global}$ represents the intercept on the log proportion scale, so we centre its prior on the typical log proportion across all locations:

$$ \beta_{global} \sim Normal(-5.5, 1) $$