##########################################################################
# Script for Rt estimation from Osiris hospital data
# - from case data (symptom onset, admission date, admission reported, infected abroad)
# - extract epicurve (admissions adjusted with nowcast for expected admissions)
# - calculate instantaneous R and case R (with 95% confidence interval)
# - save epicurve
# - save plot admission curve, epicurve and case R for Jaap
# - export json file for dashboard
##########################################################################

library(jsonlite)
library(tidyverse)
library(lubridate)
library(cowplot)
library(odbc)
library(runjags)
library(coda)

# Source functions
list.files(path = "./___UITZONDERINGSGROND_6___functions", full.names = TRUE) %>% walk(.f = source)


analysis_date <- as.Date("2020-11-20")
analysis_date <- today()

# Generate / retrieve case data
case_data <- generate_case_data_reports(date = analysis_date)

#case_data <- case_data %>% filter(no_symptoms == 0)

# Extract epicurve: expected admissions by nowcast, distinction between those infected abroad or not, and between known and unknown symptom onset
epicurve <- extract_epicurve(data = case_data, date = analysis_date, start_date = analysis_date - 100)

# distribution symptom onset to admission
SO2adm <- determine_SO2admission(data = case_data, date = analysis_date, start_date = analysis_date - 30)

plot(SO2adm, type = "l")
lines(SO2adm, col = 2)

# calculate reproduction numbers (instantaneous R and case R)
tmp <- calculate_R(epicurve = epicurve, SO2adm = SO2adm, window = 7)
epicurve <- tmp$epicurve
incidence <- tmp$incidence
saveRDS(object = epicurve, file = paste0("/___UITZONDERINGSGROND_6___epicurve_reports_", as.Date(analysis_date), ".rds"))
#saveRDS(object = incidence, file = paste0("/___UITZONDERINGSGROND_6___incidence_reports_", as.Date(analysis_date), ".rds"))

#epicurve <- readRDS(file = paste0("/___UITZONDERINGSGROND_6___epicurve_reports_", as.Date(analysis_date), ".rds"))

plot_R_Osiris(epicurve, R2plot = "tau")

# plot epicurve and reproduction number
plot_grid(plot_totcurve_Osiris(epicurve = epicurve),
          plot_Reff_Osiris(epicurve %>% mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR)), caseR = TRUE) + labs(subtitle = NULL),
          ncol = 1,
          align = "hv")

# png aan ___UITZONDERINGSGROND_2___ mailen; hij maakt presentatie voor OMT en/of Jaap

# # json file for dashboard
# reproduction.json <- epicurve %>% 
#   select(dates, caseRlower, caseR, caseRupper) %>% 
#   filter(dates >= as.Date("2020-02-17")) %>%
#   mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR),
#          caseR = round(caseR, digits = 2),
#          caseRlower = round(caseRlower, digits = 2),
#          caseRupper = round(caseRupper, digits = 2)) %>%
#   rename(Date = dates,
#          Rt_low = caseRlower,
#          Rt_avg = caseR,
#          Rt_up = caseRupper) %>%
#   toJSON 
# 
# write(x = reproduction.json, file = paste0("/___UITZONDERINGSGROND_6___COVID-19_reproductiegetal.json"))
# write(x = reproduction.json, file = paste0("/___UITZONDERINGSGROND_6___COVID-19_reproductiegetal_", as.Date(analysis_date), ".json"))

# eerste file (COVID-19_Reproductiegetal.json) op maandagmiddag (na 15:00) of dinsdagochtend (voor 13:00) op R:\Projecten\COVID-19\Surveillance\Open_data zetten


epicurve %>% filter(dates == analysis_date - 14) %>% select(dates, caseR, caseRlower, caseRupper)

ggplot(data = epicurve %>% filter(dates >= as.Date("2020-07-01")) %>% mutate(fracabroad = admission_infectedabroad/admission),
       mapping = aes(x = dates, y = fracabroad)) + 
  geom_line() + 
  theme_minimal()

# plot epicurve and reproduction number
plot_grid(plot_totcurve_Osiris_reports(epicurve),
          plot_Reff_Osiris(epicurve %>% mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR)), caseR = TRUE) +
            labs(subtitle = paste("OSIRIS data", last(epicurve$dates))) +
            coord_cartesian(
              ylim = c(0,2)) +
            scale_x_date(
              limits = c(as.Date("2020-06-12"), last(epicurve$dates)+1),
              expand = expansion(add = -0.5),
              breaks = "4 days",
              minor_breaks = "1 day",
              date_labels = "%e %b") +
            labs(subtitle = NULL),
          ncol = 1,
          align = "hv") 

ggsave(filename = paste0("/___UITZONDERINGSGROND_6___epicurveReff_Osiris_reports_", as.Date(analysis_date), ".png"), width = 7, height = 6, dpi = 300)


ggplot(data = case_data %>% filter(admission_reported >= as.Date("2020-07-01") & !is.na(admission_date)) %>% 
         mutate(delay = as.integer(admission_reported - admission_date)) %>% 
         group_by(admission_reported) %>% 
         summarise(delay_mean = mean(delay),
                   delay_lower = quantile(delay, 0.025),
                   delay_upper = quantile(delay, 0.975)),
       mapping = aes(x = admission_reported, y = delay_mean, ymin = delay_lower, ymax = delay_upper)) +
  geom_ribbon(col = NA, fill = adjustcolor(4, alpha = 0.3)) +
  geom_point(col = 4) +
  labs(x = "datum fiattering melding",
       y = "aantal dagen tussen start en fiattering melding") +
  coord_cartesian(ylim = c(0, 7)) +
  theme_minimal()

ggsave(filename = paste0("/___UITZONDERINGSGROND_6___delay_Osiris_reports_", as.Date(analysis_date), ".png"), width = 7, height = 6, dpi = 300)

ggplot(data = case_data %>% filter(!is.na(symptom_onset)) %>% 
         mutate(delay = as.integer(admission_date - symptom_onset)) %>% 
         group_by(admission_date) %>% 
         summarise(delay_mean = mean(delay),
                   delay_lower = quantile(delay, 0.025),
                   delay_upper = quantile(delay, 0.975)),
       mapping = aes(x = admission_date, y = delay_mean, ymin = delay_lower, ymax = delay_upper)) +
  geom_ribbon(col = NA, fill = adjustcolor("#01689b", alpha = 0.3)) +
  geom_point(col = "#01689b") +
  labs(x = "datum melding",
       y = "aantal dagen tussen EZD en melding") +
  coord_cartesian(xlim = c(as.Date("2020-07-01"), NA), ylim = c(0, 14)) +
  theme_minimal()


epicurve %>% 
  filter(dates > as.Date("2020-06-12")) %>% 
  select(dates, admission, incidence_mean, caseR) %>% 
  rename(reported = admission, symptom_onset = incidence_mean) %>% 
  writexl::write_xlsx("./___UITZONDERINGSGROND_6___epicurve_Jaap_2020-11-20.xlsx")


# omit asymptomatic cases

files <- list.files("/___UITZONDERINGSGROND_6___Previous", full.names = TRUE)
files <- c(files, list.files("/___UITZONDERINGSGROND_6___Geschoond", full.names = TRUE))
files <- files[grepl(files, pattern = "rds")]
files <- tibble(file_name = files, 
                date = files %>% str_extract_all(pattern = "(Data_[0-9]+)") %>% unlist %>% gsub(pattern = "Data_", replacement = "") %>% as.Date(format = "%Y%m%d"),
                time = files %>% file.info %>% pull(ctime))

file_name <- files %>% 
    filter(date == analysis_date) %>%
    filter(time == max(time)) %>% # get most recent file (if multiple files exist) 
    pull(file_name)
  
dataOsiris <- readRDS(file_name)

case_data <- dataOsiris %>% 
  filter(., !(EIGENAARDesc %in% c("Sint Eustatius", "Bonaire", "Saba", "Sint Maarten", "Aruba", "Curaçao", "Curacao"))) %>% 
  filter(NCOVVast1eziektedag != "NVT") %>% select(OSIRISNR) %>% inner_join(case_data)