In this document, we assess our predictions of a simple definition of direction. First we get the data.

# Load the evaluations
save_data <- readRDS(file.path(analysis_dir, paste(analysis_type,"prediction.rds",sep = "-")))
evals <- save_data$evals
evals <- evals %>% rename(target_date = target_end_date)
forecasting_task <- save_data$forecasting_task
forecaster_args <- save_data$forecaster_args

# Load the trained models
forecasters <- unique(evals$forecaster)
trained_models <- lapply(forecasters, FUN = function(forecaster){
  readRDS_from_dir(file.path(model_dir,forecaster))
})
names(trained_models) <- forecasters

# Create coefficient dfs
coefs <- lapply(trained_models,create_coef_df) %>% 
  bind_rows(.id = "forecaster") %>%
  mutate(feature_name = parse_feature_name(feature_name))

# Omit whatever was asked to be omitted.
evals <- evals %>% filter(!(forecaster %in% omitted_forecasters) & 
                          !(ahead %in% omitted_aheads) &
                          !(forecast_date %in% omitted_forecast_dates))
coefs <- coefs %>% filter(!(forecaster %in% omitted_forecasters) & 
                      !(ahead %in% omitted_aheads) &
                      !(forecast_date %in% omitted_forecast_dates))

Our forecasting task is:

Geographical resolution: state
Response variable: hhs, confirmed_admissions_covid_1d_7dav
Incidence period: day
Forecast date range: 2020-10-01 to 2021-11-01
Number of days ahead forecast: 7 to 28.

We transform the response into a binary variable for classification. We will call this variable direction. The parameters used to define direction are

upswing threshold (ratio of present - ast value over past value: .1
minimum threshold (minimum number of cases needed at past value in order to qualify) 130/7.

Let \(N_i^{t}\) denote the response variable, observed at a given {location = i,day = t}. A given {location = i,forecast_date = t,ahead = a} triplet qualifies as an increase if

\[ \frac{{N}_i^{t + a} - {N}_i^{t}}{{N}_{i}^{t}} \geq s_0~~\textrm{and}~~ \bar{N}_i^{t} \geq t_0 \] where \(s_0\) is the upswing threshold, and \(t_0\) is the minimum threshold. If \(N_i^{t} < t_0\), then the response is coded as NA (i.e. thrown out).

Data wrangling

We restrict outselves to evaluation on a “reasonable” subset of {location,forecast_date,forecaster} triplets.

## Handle missing predictions.

# Filter to common forecast dates
common_forecast_dates <- find_common_forecast_dates(evals)
evals <- evals %>% filter(forecast_date %in% common_forecast_dates)
coefs <- coefs %>% filter(forecast_date %in% common_forecast_dates)

# Now, we do the follow operations in order.
# 
# (1) We get rid of any forecaster which is NA 
#   for more than max_forecaster_na of all tasks.
#
# (2) We get rid of any location for which a remaining forecaster has been NA
#   for more than max_location_na of all tasks.
#
# (3) We get rid of any remaining tasks at which a remaining forecaster has been
#   NA.
max_forecaster_na <- 1
max_location_na <- 1
n_forecast_date <- length(unique(evals$forecast_date))
n_ahead <- length(unique(evals$ahead))
n_locations <- length(unique(evals$geo_value))

na_tasks <- evals %>% 
  filter(is.na(value) | is.na(actual)) %>%
  select(geo_value, forecaster, forecast_date, ahead)

na_prop_by_forecaster <- na_tasks %>%
  group_by(forecaster) %>% 
  summarize(prop = n()/(n_forecast_date * n_ahead * n_locations), .groups = "drop")
na_forecasters <- na_prop_by_forecaster %>%
  filter(prop > max_forecaster_na) %>% 
  pull(forecaster) %>% unique()
na_prop_by_location <- na_tasks %>%
  group_by(forecaster,geo_value) %>%
  summarize(prop = n() / (n_ahead * n_forecast_date),.groups = "drop") %>%
  filter(!(forecaster %in% na_forecasters))
na_locations <- na_prop_by_location %>%  
  filter(prop > max_location_na) %>% 
  pull(geo_value) %>% unique()
na_tasks <- na_tasks %>%
  filter(!(forecaster %in% na_forecasters | geo_value %in% na_locations ))%>%
  select(-forecaster) %>% distinct()

evals <- evals %>%
  filter(!(forecaster %in% na_forecasters | geo_value %in% na_locations)) %>%
  anti_join(na_tasks)

Now we aggregate our evaluations (using power at different levels of type I error, and area under the ROC curve), and our estimated coefficients (using median). Aggregations are over all time and per forecast yearmonth.

### Aggregations over all forecast dates.

# ROC
roc_by_alpha <- function(x,alpha = seq(.01,1,by = .02)){
  power <- roc(x$value,x$actual,alpha)
}

all_day_roc <- evals %>%
  select(forecaster,ahead,forecast_date,value,actual) %>%
  group_by(forecaster,ahead) %>%
  group_modify( ~ roc_by_alpha(.x))

# AUC
all_day_evals <- evals %>% 
  select(forecaster,ahead,forecast_date,value,actual) %>%
  group_by(forecaster,ahead) %>%
  summarize(auc = auc(value,actual,alpha = seq(0,1,by = .01)),.groups = "drop")

# Coefficients
all_day_coefs <- coefs %>%
  group_by(ahead, feature_name, forecaster) %>% 
  summarise(value = median(value), .groups = "drop")

### Aggregation per forecast-year month.

# ROC
per_yearmonth_roc <- evals %>%
  select(forecaster,ahead,forecast_date,value,actual) %>%
  mutate(forecast_month = lubridate::month(forecast_date),
         forecast_year = lubridate::year(forecast_date)) %>%
  group_by(forecaster,ahead,forecast_month,forecast_year) %>%
  group_modify( ~ roc_by_alpha(.x))

# Compute AUC over each month
per_yearmonth_evals <- evals %>% 
  select(forecaster,ahead,forecast_date,value,actual) %>%
  mutate(forecast_month = lubridate::month(forecast_date),
         forecast_year = lubridate::year(forecast_date)) %>%
  group_by(forecaster,ahead,forecast_month,forecast_year) %>%
  summarize(auc = auc(value,actual,alpha = seq(0,1,by = .01)),.groups = "drop")

# Coefficients
per_yearmonth_coefs <- coefs %>%
  mutate(forecast_month = lubridate::month(forecast_date),
         forecast_year = lubridate::year(forecast_date)) %>%
  group_by(ahead, feature_name, forecaster,forecast_month,forecast_date) %>% 
  summarise(value = median(value), .groups = "drop")

EDA

The proportion of places that qualify as an increase in direction, over time.

# Plot the proportion of increases over time.
df <- evals %>% 
  filter(forecaster == "AR3") %>%
  rename(direction = actual) %>%
  select(geo_value,target_date,direction) %>%
  group_by(target_date) %>%
  summarize(prop = sum(direction,na.rm = TRUE)/n(), .groups = "drop")

ggplot(df, aes(x = target_date, y = prop)) +
  geom_line() +
  geom_point() + 
  scale_x_date(date_minor_breaks = "1 month") +
  labs(subtitle = subtitle)

Evaluation

First up, ROC curves by ahead, aggregated over all days.

ROC curves

7 days ahead

ggplot(all_day_roc %>% filter(ahead == 7), 
       aes(x = alpha, y = result, color = forecaster)) +
  geom_point() + 
  geom_line() +
  labs(subtitle = subtitle,
       x = "False positive rate",
       y = "True positive rate")

10 days ahead

ggplot(all_day_roc %>% filter(ahead == 10), 
       aes(x = alpha, y = result, color = forecaster)) +
  geom_point() + 
  geom_line() +
  labs(subtitle = subtitle,
       x = "False positive rate",
       y = "True positive rate")

14 days ahead

ggplot(all_day_roc %>% filter(ahead == 14), 
       aes(x = alpha, y = result, color = forecaster)) +
  geom_point() + 
  geom_line() +
  labs(subtitle = subtitle,
       x = "False positive rate",
       y = "True positive rate")

21 days ahead

ggplot(all_day_roc %>% filter(ahead == 21), 
       aes(x = alpha, y = result, color = forecaster)) +
  geom_point() + 
  geom_line() +
  labs(subtitle = subtitle,
       x = "False positive rate",
       y = "True positive rate")

Next up, AUC by ahead, aggregated over all dates. The tabs differ in the notion of aggregation:

In the first tab, we aggregate by taking AUC over all forecast dates.
In the second tab, we group by forecast year and month, then take AUC over all forecast dates.

AUC curves by ahead

All time

ggplot(all_day_evals, aes(x = ahead, y = auc, color = forecaster)) +
  geom_point() +
  geom_line() +
  labs(subtitle = subtitle)

Per Month

ggplot(per_yearmonth_evals, aes(x = ahead, y = auc, color = forecaster)) +
  geom_point() +
  geom_line() +
  labs(subtitle = subtitle) +
  facet_wrap(vars(forecast_year,forecast_month), scales = "free")

Finally, difference in fitted values between forecasters, depending on the truth.

Fitted values

All aheads

ggplot(diff_evals, aes(x = diff, fill = forecaster)) +
  geom_histogram(bins = 200) +
  labs(subtitle = subtitle) +
  facet_grid(rows = vars(actual),
             cols = vars(forecaster),
             scales = "fixed") +
  geom_vline(xintercept = 0)

By ahead

7 days ahead

ggplot(diff_evals %>% filter(ahead == 7), aes(x = diff, fill = forecaster)) +
  geom_histogram(bins = 200) +
  labs(subtitle = subtitle) +
  facet_grid(rows = vars(actual),
             cols = vars(forecaster),
             scales = "free_x") +
  geom_vline(xintercept = 0)

14 days ahead

ggplot(diff_evals %>% filter(ahead == 14), aes(x = diff, fill = forecaster)) +
  geom_histogram(bins = 200) +
  labs(subtitle = subtitle) +
  facet_grid(rows = vars(actual),
             cols = vars(forecaster),
             scales = "free_x") +
  geom_vline(xintercept = 0)

21 days ahead

ggplot(diff_evals %>% filter(ahead == 21), aes(x = diff, fill = forecaster)) +
  geom_histogram(bins = 200) +
  labs(subtitle = subtitle) +
  facet_grid(rows = vars(actual),
             cols = vars(forecaster),
             scales = "free_x") +
  geom_vline(xintercept = 0)

Estimates

Estimated coefficients by value of ahead, aggregated over all dates.

ggplot(all_day_coefs, aes(x = ahead, y = value,color = forecaster)) + 
  geom_point() + 
  geom_line() +
  facet_wrap(vars(feature_name), scales = "free", ncol = 3)

Now, estimated coefficients over time, for ahead = c(7,10,14,21).

Parameter estimates over time

Unsmoothed

7 days ahead

ggplot(coefs %>% filter(ahead == 7), aes(x = forecast_date, y = value, color = forecaster)) + 
    geom_point(cex = .25) +
    geom_line() +
    scale_x_date(date_minor_breaks = "1 month") +
    facet_wrap(vars(feature_name), scales = "free", ncol = 3)

10 days ahead

ggplot(coefs %>% filter(ahead == 10), aes(x = forecast_date, y = value, color = forecaster)) + 
    geom_point(cex = .25) +
    geom_line() +
    scale_x_date(date_minor_breaks = "1 month") + 
    facet_wrap(vars(feature_name), scales = "free", ncol = 3)

14 days ahead

ggplot(coefs %>% filter(ahead == 14), aes(x = forecast_date, y = value, color = forecaster)) + 
    geom_point(cex = .25) +
    geom_line() +
    scale_x_date(date_minor_breaks = "1 month") +
    facet_wrap(vars(feature_name), scales = "free", ncol = 3)

21 days ahead

ggplot(coefs %>% filter(ahead == 21), aes(x = forecast_date, y = value, color = forecaster)) + 
    geom_point(cex = .25) +
    geom_line() +
    scale_x_date(date_minor_breaks = "1 month") +
    facet_wrap(vars(feature_name), scales = "free", ncol = 3)

1. Produce and evaluate direction predictions for COVID hospitalizations

Data wrangling

EDA

Evaluation

ROC curves

7 days ahead

10 days ahead

14 days ahead

21 days ahead

AUC curves by ahead

All time

Per Month

Fitted values

All aheads

By ahead

7 days ahead

14 days ahead

21 days ahead

Estimates

Parameter estimates over time

Unsmoothed

7 days ahead

10 days ahead

14 days ahead

21 days ahead