```
<- readRDS(url(
us_raw_history_dfs "https://github.com/cmu-delphi/epipredict/raw/dev/vignettes/articles/all_states_covidcast_signals.rds"
))
<- us_raw_history_dfs[[1]] %>%
us_cli_archive select(geo_value, time_value, version = issue, percent_cli = value) %>%
as_epi_archive(compactify = TRUE)
<- us_raw_history_dfs[[2]] %>%
us_cases_archive select(geo_value, time_value, version = issue, case_rate = value) %>%
as_epi_archive(compactify = TRUE)
<- epix_merge(
us_archive
us_cli_archive, us_cases_archive,sync = "locf", compactify = TRUE
)
```

# 14 Sliding version-unaware and version-aware ARX forecasters across dates

A key function from the epiprocess package is `epix_slide()`

(refer to the following vignette for the basics of the function: “Work with archive objects and data revisions”) which allows performing version-aware computations. That is, the function only uses data that would have been available as of time t for that reference time.

In this vignette, we use `epix_slide()`

for backtesting our `arx_forecaster`

on historical COVID-19 case data from the US and from Canada. We first examine the results from a version-unaware forecaster, comparing two different fitting engines and then we contrast this with version-aware forecasting. The former will proceed by constructing an `epi_archive`

that erases its version information and then use `epix_slide()`

to forecast the future. The latter will keep the versioned data and proceed similarly by using `epix_slide()`

to forecast the future.

## 14.1 Version-unaware forecasting

### 14.1.1 Example using CLI and case data from US states

First, we download the version history (i.e. archive) of the percentage of doctor’s visits with CLI (COVID-like illness) computed from medical insurance claims and the number of new confirmed COVID-19 cases per 100,000 population (daily) for all 50 states from the COVIDcast API. We process as before, with the modification that we use `sync = "locf"`

in `epix_merge()`

so that the last version of each observation can be carried forward to extrapolate unavailable versions for the less up-to-date input archive.

We then get latest snapshot of the data from the archive by using `epix_as_of()`

. We then create fake version information by setting `version = time_value`

. This creates an archive that pretends to have the latest data available (since at version time `x`

it has all the data up to time_value `x`

, which in reality is unrealistic because the time values of the data received at version time `x`

often lags by a few days, not to mention the later corrections that are amended to the data).

```
# Get latest snapshot of data and pretend it's an archive
<- us_archive %>%
us_latest epix_as_of(version = max(.$versions_end)) %>%
mutate(version = time_value) %>%
as_epi_archive()
<- seq(
fc_time_values from = as.Date("2020-08-01"),
to = as.Date("2021-11-01"),
by = "1 month"
)<- c(7, 14, 21, 28)
aheads
<- function(epi_archive, outcome, predictors, ahead = 7, engine) {
forecast_k_week_ahead %>%
epi_archive epix_slide(
.f = function(x, gk, rtv) {
arx_forecaster(
x, outcome, predictors, engine,args_list = arx_args_list(ahead = ahead)
$predictions %>%
)mutate(engine_type = engine$engine) %>%
pivot_quantiles_wider(.pred_distn)
},.before = 120,
.versions = fc_time_values
)
}
# Generate the forecasts and bind them together
<- bind_rows(
forecasts_version_unaware map(aheads, ~ forecast_k_week_ahead(
us_latest,outcome = "case_rate",
predictors = c("case_rate", "percent_cli"),
ahead = .x,
engine = linear_reg()
)),map(aheads, ~ forecast_k_week_ahead(
us_latest,outcome = "case_rate",
predictors = c("case_rate", "percent_cli"),
ahead = .x,
engine = rand_forest(mode = "regression")
)) )
```

Here, `arx_forecaster()`

does all the heavy lifting. It creates leads of the target (respecting time stamps and locations) along with lags of the features (here, the response and doctors visits), estimates a forecasting model using the specified engine, creates predictions, and non-parametric confidence bands.

To see how the predictions compare, we plot them on top of the latest case rates. Note that even though we’ve fitted the model on all states, we’ll just display the results for two states, California (CA) and Florida (FL), to get a sense of the model performance while keeping the graphic simple.

## Code

```
<- forecasts_version_unaware %>%
forecasts_filtered tibble() %>%
filter(geo_value %in% c("ca", "fl"))
<- us_latest$DT %>%
latest_data_filtered tibble() %>%
filter(geo_value %in% c("ca", "fl"))
ggplot(forecasts_filtered, aes(x = target_date, group = forecast_date, fill = engine_type)) +
geom_line(
data = latest_data_filtered, aes(x = time_value, y = case_rate),
inherit.aes = FALSE, color = "gray50"
+
) geom_ribbon(aes(ymin = `0.05`, ymax = `0.95`), alpha = 0.4) +
geom_line(aes(y = .pred)) +
geom_point(aes(y = .pred), size = 0.5) +
geom_vline(aes(xintercept = forecast_date), linetype = 2, alpha = 0.5) +
facet_grid(engine_type ~ geo_value, scales = "free") +
scale_x_date(minor_breaks = "month", date_labels = "%b %y") +
scale_fill_brewer(palette = "Set1") +
scale_y_continuous(expand = expansion(c(0, 0.05))) +
labs(x = "Date", y = "Reported COVID-19 case rates") +
theme(legend.position = "none")
```

For the two states of interest, simple linear regression clearly performs better than random forest in terms of accuracy of the predictions and does not result in such in overconfident predictions (overly narrow confidence bands). Though, in general, neither approach produces amazingly accurate forecasts. This could be because the behaviour is rather different across states and the effects of other notable factors such as age and public health measures may be important to account for in such forecasting. Including such factors as well as making enhancements such as correcting for outliers are some improvements one could make to this simple model.^{1}

### 14.1.2 Example using case data from Canada

By leveraging the flexibility of `epiprocess`

, we can apply the same techniques to data from other sources. Since some collaborators are in British Columbia, Canada, we’ll do essentially the same thing for Canada as we did above.

The COVID-19 Canada Open Data Working Group collects daily time series data on COVID-19 cases, deaths, recoveries, testing and vaccinations at the health region and province levels. Data are collected from publicly available sources such as government datasets and news releases. Unfortunately, there is no simple versioned source, so we have created our own from the Github commit history.

First, we load versioned case rates at the provincial level. After converting these to 7-day averages (due to highly variable provincial reporting mismatches), we then convert the data to an `epi_archive`

object, and extract the latest version from it. Finally, we run the same forcasting exercise as for the American data, but here we compare the forecasts produced from using simple linear regression with those from using boosted regression trees.

```
# source("drafts/canada-case-rates.R)
<- epidatasets::can_prov_cases
can <- can %>%
can group_by(version, geo_value) %>%
arrange(time_value) %>%
mutate(cr_7dav = RcppRoll::roll_meanr(case_rate, n = 7L)) %>%
as_epi_archive(compactify = TRUE)
<- epix_as_of(can, max_version = max(can$DT$version)) %>%
can_latest mutate(version = time_value) %>%
as_epi_archive()
# Generate the forecasts, and bind them together
<- bind_rows(
can_fc map(
aheads,~ forecast_k_week_ahead(can_latest, "cr_7dav", "cr_7dav", .x, linear_reg())
),map(
aheads,~ forecast_k_week_ahead(can_latest, "cr_7dav", "cr_7dav", .x, boost_tree(mode = "regression", trees = 20))
) )
```

The first figure shows the results for all of the provinces using linear regression.

## Code

```
ggplot(
%>% filter(engine_type == "lm"),
can_fc aes(x = target_date, group = forecast_date)
+
) coord_cartesian(xlim = lubridate::ymd(c("2020-12-01", NA))) +
geom_line(
data = can_latest$DT %>% tibble(), aes(x = time_value, y = cr_7dav),
inherit.aes = FALSE, color = "gray50"
+
) geom_ribbon(aes(ymin = `0.05`, ymax = `0.95`, fill = geo_value),
alpha = 0.4
+
) geom_line(aes(y = .pred)) +
geom_point(aes(y = .pred), size = 0.5) +
geom_vline(aes(xintercept = forecast_date), linetype = 2, alpha = 0.5) +
facet_wrap(~geo_value, scales = "free_y", ncol = 3) +
scale_x_date(minor_breaks = "month", date_labels = "%b %y") +
scale_y_continuous(expand = expansion(c(0, 0.05))) +
labs(
title = "Using simple linear regression", x = "Date",
y = "Reported COVID-19 case rates"
+
) theme(legend.position = "none")
```

Compare those forecasts with a related set using Gradient Boosting.

## Code

```
ggplot(
%>% filter(engine_type == "xgboost"),
can_fc aes(x = target_date, group = forecast_date)
+
) coord_cartesian(xlim = lubridate::ymd(c("2020-12-01", NA))) +
geom_line(
data = can_latest$DT %>% tibble(), aes(x = time_value, y = cr_7dav),
inherit.aes = FALSE, color = "gray50"
+
) geom_ribbon(aes(ymin = `0.05`, ymax = `0.95`, fill = geo_value),
alpha = 0.4
+
) geom_line(aes(y = .pred)) +
geom_point(aes(y = .pred), size = 0.5) +
geom_vline(aes(xintercept = forecast_date), linetype = 2, alpha = 0.5) +
facet_wrap(~geo_value, scales = "free_y", ncol = 3) +
scale_x_date(minor_breaks = "month", date_labels = "%b %y") +
scale_y_continuous(expand = expansion(c(0, 0.05))) +
labs(
title = "Using boosted regression trees", x = "Date",
y = "Reported COVID-19 case rates"
+
) theme(legend.position = "none")
```

Both approaches tend to produce quite volatile forecasts (point predictions) and/or are overly confident (very narrow bands), particularly when boosted regression trees are used. But as this is meant to be a simple demonstration of sliding with different engines in `arx_forecaster`

, we may devote another vignette to work on improving the predictive modelling using the suite of tools available in epipredict.

## 14.2 Version-aware forecasting

### 14.2.1 Example using case data from US states

We will now run pseudoprospective forecasts based on properly-versioned data (that would have been available in real-time) to forecast future COVID-19 case rates from current and past COVID-19 case rates for all states. All we have to do is use the historical archive of the data with version information, `us_archive`

, instead of `us_latest`

like we did above, in the argument to our forecaster wrapper `forecast_k_week_ahead()`

. Below we do that computation, tag it, and combine it with the forecasts from one of the engines made above.

```
<- map(aheads, ~ forecast_k_week_ahead(
forecasts_version_aware
us_archive,outcome = "case_rate",
predictors = c("case_rate", "percent_cli"),
ahead = .x,
engine = linear_reg()
%>%
)) bind_rows() %>%
mutate(version = "version faithful")
```

Now we can plot the results on top of the latest case rates. As before, we will only display and focus on the results for FL and CA for simplicity.

## Code

```
<- bind_rows(
forecasts_filtered
forecasts_version_aware,%>%
forecasts_version_unaware filter(engine_type == "lm") %>%
mutate(version = "version unfaithful")
%>%
) tibble() %>%
filter(geo_value %in% c("ca", "fl"))
<- us_latest$DT %>%
latest_data_filtered tibble() %>%
select(-version) %>%
filter(geo_value %in% c("ca", "fl"))
ggplot(forecasts_filtered, aes(x = target_date, group = forecast_date, fill = version)) +
geom_line(
data = latest_data_filtered, aes(x = time_value, y = case_rate),
inherit.aes = FALSE, color = "gray50"
+
) geom_ribbon(aes(ymin = `0.05`, ymax = `0.95`), alpha = 0.4) +
geom_line(aes(y = .pred)) +
geom_point(aes(y = .pred), size = 0.5) +
geom_vline(aes(xintercept = forecast_date), linetype = 2, alpha = 0.5) +
facet_grid(version ~ geo_value, scales = "free") +
scale_x_date(minor_breaks = "month", date_labels = "%b %y") +
scale_fill_brewer(palette = "Set1") +
scale_y_continuous(expand = expansion(c(0, 0.05))) +
labs(x = "Date", y = "Reported COVID-19 case rates") +
theme(legend.position = "none")
```

Again, we observe that the results are not great for these two states, but that’s likely due to the simplicity of the model (ex. the omission of key factors such as age and public health measures) and the quality of the data (ex. we have not personally corrected for anomalies in the data).

We shall leave it to the reader to try the above version aware and unaware forecasting exercise on the Canadian case rate data. The above code for the American state data should be readily adaptable for this purpose.

Note that, despite the above caveats, simple models like this tend to out-perform many far more complicated models in the online Covid forecasting due to those models high variance predictions.↩︎