Script and Analytical Data Set for forecasting Covid in Denmark

forecast_ads.csv

script.R

library(fable)
library(tsibble)
library(tsibbledata)
library(lubridate)
library(fable.prophet)
library(tidyverse)

################################################################################

# Read in ADS data from csv gist (forecast_ads.csv)
# ads <- read_csv(PATH)

fit <- ads %>%
  filter(date >= "2020-08-03") %>%
  stretch_tsibble(.step = 7, .init = 70) %>%
  model(
    mean = MEAN(rm),
    ets = ETS(rm),
    arima = ARIMA(rm),
    naive = NAIVE(rm),
    nnetar = NNETAR(rm),
    tslm = TSLM(rm),
    prophet = prophet(rm)
  ) 

fit %>%
  glance()

f <- fit %>%
  forecast(h = "2 weeks")

f

ac <- accuracy(f, ads) %>%
  arrange(RMSE)

ac %>%
  arrange(regionnavn, MAPE) %>%
  group_by(regionnavn) %>%
  slice(1:3) %>%
  ungroup() %>%
  pull(.model) %>%
  unique()

saveRDS(f, "f.RDS")

################################################################################

weekstarts <- ymd("2020-10-26") + c(0, 7, 14, 21, 28)

weekstart <- weekstarts[1]

out <- map_dfr(weekstarts, function(weekstart){
  model_data <- ads %>%
    filter(date >= "2020-08-03",
           date < weekstart)
  
  fit2 <-  model_data %>%
    model(
      naive = NAIVE(rm),
      ets = ETS(rm),
      nnetar = NNETAR(rm),
      arima = ARIMA(rm)
    ) 
  
  f2 <- fit2 %>%
    forecast(h = "2 weeks")
  
  f2 %>%
    as_tibble() %>%
    mutate(id = str_c(weekstart, "_", .model))
})

saveRDS(out, "forecasts.RDS")

################################################################################

out <- readRDS("~/forecast/forecasts.RDS")
f <- readRDS("~/forecast/f.RDS")

ac <- accuracy(f, ads) %>%
  arrange(RMSE)

weights <- ac %>%
  arrange(regionnavn, MAPE) %>%
  select(.model, regionnavn, MAPE) 

en <- out %>%
  select(-rm) %>%
  arrange(regionnavn, date) %>%
  left_join(weights) %>%
  mutate(weekstart = id %>% str_sub(1,10)) %>%
  group_by(regionnavn, date, weekstart) %>%
  mutate(inv = 1/MAPE,
         w = inv / sum(inv)) %>%
  summarise(.mean = weighted.mean(.mean, w)) %>%
  ungroup() %>%
  mutate(.model = "ensemble", 
         id = str_c(weekstart, "_", .model)) %>%
  select(-weekstart)

forecastdata <- out %>% select(regionnavn, .model, date, .mean, id)

mfc <- forecastdata %>%
  bind_rows(en)

region_ads <- ads %>% filter(regionnavn == "Hovedstaden")
region_fc <- mfc %>% filter(regionnavn == "Hovedstaden")

N = 20000

ggplot() + 
  geom_line(data = region_ads %>% filter(date >= "2020-10-01"), 
            aes(date, pospct * N)) + 
  geom_line(data = region_fc %>% 
              filter(.model == "ensemble") %>%
              filter(id == max(id)), 
            aes(date, .mean * N, color = .model, group = id)) + 
  theme_minimal()

####

test_df <- tribble(
  ~regionnavn, ~antal,
  "Hovedstaden", 20000,
  "Midtjylland", 10000,
  "Nordjylland", 6000,
  "Sjælland", 6000,
  "Syddanmark", 9000
)

facet_labs = glue::glue("{test_df$regionnavn} med {test_df$antal} test") %>%
  as.character()
names(facet_labs) = test_df$regionnavn

myfc = mfc %>%
  left_join(test_df) %>%
  mutate(cases = .mean * antal) %>% 
  filter(.model == "ensemble") %>%
  filter(id == max(id)) %>%
  select(date, regionnavn, cases)

myads = ads %>% 
  left_join(test_df) %>%
  filter(date >= "2020-08-01") %>%
  mutate(cases = pospct * antal,
         cases_rm = rm * antal) %>%
  select(date, regionnavn, cases, cases_rm)

ggplot() + 
  geom_point(data = myads, aes(date, cases, group = regionnavn),
             shape = 1) + 
  geom_line(data = myads, aes(date, cases_rm, group = regionnavn),
            color = "red") + 
  geom_line(data = myfc, aes(date, cases, group = regionnavn), 
            color = "red", linetype = "dashed") + 
  theme_minimal() + 
  facet_wrap(~regionnavn, scales = "free_y", ncol = 2,
             labeller = labeller(regionnavn = facet_labs))

ggsave("regional_forecas.png", width = 10, height = 10)


save(myfc, myads, file = "uge48_forecast.Rdata")