###########################################################
# Epicurve symptom onset of hospitalised cases in Osiris data
# 
# 20-03-2020 ___UITZONDERINGSGROND_2___
###########################################################

library(tidyverse)
library(lubridate)
library(odbc)
library(runjags)
library(coda)
library(cowplot)
library(gridGraphics)
#ibrary(igraph)
#library(shiny)
#library(visNetwork)
#library(readxl)
#library(writexl)

list.files(path = "Functies", full.names = TRUE) %>% 
  walk(.f = source)


# get most recent case data from Osiris
files <- list.files("/___UITZONDERINGSGROND_6___Data_OSIRIS_geschoond", full.names = TRUE)
files <- files[grepl(files, pattern = "rds")]
file.name <- files[files %>% file.info %>% pull(ctime) %>% which.max]
dataOsiris <- read_rds(file.name)


data <- dataOsiris %>% 
  #filter(Leeftijd >= 70) %>%
  mutate(CREATIEDT     = CREATIEDT %>% as.Date(format = "%d-%m-%Y"),
         NCOVdat1ezkhopn = NCOVdat1ezkhopn %>% as.Date(format = "%Y-%m-&d")) %>% 
  filter(NCOVpatZhs == "J")

start.date <- as.Date("2020-02-15")
date.seq <- seq(from = start.date, to = today(), by = "day")

# Run the nowcast
nowcast_data <- nowcast(
  case.data = data %>%
    transmute(onset.date = ZIE1eZiekteDt, report.date = CREATIEDT) %>% 
    drop_na,
  start.date   = start.date,
  nowcast.date = today())


# # Run the nowcast
# nowcast_data <- nowcast(
#   case.data = data %>%
#     transmute(onset.date = NCOVdat1ezkhopn, report.date = CREATIEDT) %>% 
#     drop_na,
#   start.date   = start.date,
#   nowcast.date = today())


# Assign to plot_incidenceziektedag, so it can be re-used for the nowcast plot
plot_incidenceziektedag <- ggplot(
  data = data %>%
    mutate(ZIE1eZiekteDt = ZIE1eZiekteDt %>% factor(levels = date.seq %>% as.character)) %>%
    group_by(ZIE1eZiekteDt, .drop = FALSE) %>%
    summarize(Aantal = n()) %>%
    drop_na %>%
    mutate(ZIE1eZiekteDt = ZIE1eZiekteDt %>% as.Date),
  mapping = aes(x = ZIE1eZiekteDt, y = Aantal)) +
  geom_col(fill = "#01689b", width = 0.9) +
  labs(x = "Datum eerste ziektedag", y = "Aantal") +
  theme_minimal()

# Epicurve Nowcast - Totaal aantal cases naar eerste ziektedag
plot_incidenceziektedag +
  geom_crossbar(
    data = nowcast_data,
    mapping = aes(x = onset.date, y = nowcast_p50, ymin = nowcast_p10, ymax = nowcast_p90),
    fill = adjustcolor("#e17000", alpha = 0.5),
    colour = NA,
    width = 0.9) +
  geom_crossbar(
    data = nowcast_data,
    mapping = aes(x = onset.date, y = nowcast_p50, ymin = nowcast_p25, ymax = nowcast_p75),
    fill = adjustcolor("#e17000", alpha = 0.5),
    colour = NA,
    width = 0.9) +
  # geom_ribbon(
  #   mapping = aes(ymin = trend_lwr, ymax = trend_upr),
  #   fill = "black",
  #   alpha = 0.1) +
  geom_line(
    data = nowcast_data,
    mapping = aes(x = onset.date, y = trend_med),
    colour = "black",
    size = 0.25) +
  scale_x_date(
    breaks = "1 week",
    minor_breaks = "1 day",
    date_labels = "%d %b") +
  scale_y_continuous(
    breaks = c(0, outer(c(1, 2, 5), 10^(0:3)) %>% as.vector),
    minor_breaks = outer(c(3, 4, 6, 7, 8, 9), 10^(0:3)) %>% as.vector) +
  coord_trans(y = "log1p")





# median delay SO to reporting: 6 days
interval <- as.Date(data$CREATIEDT, "%Y-%m-%d") - data$ZIE1eZiekteDt
table(interval, useNA = "ifany")
hist(as.integer(interval), breaks = -1:30)
median(interval, na.rm = TRUE)
test <- ecdf(interval)
plot(test(0:30))



# median delay SO to hospital: 6 days
interval <- as.Date(data$NCOVdat1ezkhopn, "%Y-%m-%d") - data$ZIE1eZiekteDt
interval <- interval[interval>-40]
table(interval, useNA = "ifany")
hist(as.integer(interval), breaks = -32:23)
median(interval[interval > 0], na.rm = TRUE)

# median delay hospital to reporting: 1 day
interval <- data$MELGGDOntvDt - as.Date(data$NCOVdat1ezkhopn, "%Y-%m-%d")
interval <- interval[interval < 40]
table(interval, useNA = "ifany")
hist(as.integer(interval), breaks = -5:33)
median(interval, na.rm = TRUE)


# fraction with missing SO
sum(is.na(data$ZIE1eZiekteDt))/nrow(data)

# distribution of missing SO over reporting date -> 
ggplot(data = data %>% filter(Leeftijd >= 70),
       mapping = aes(x = CREATIEDT, fill = is.na(ZIE1eZiekteDt))) +
  geom_bar(position = "fill") + 
  theme_minimal()

# epicurve of SO of hospitalised cases
ggplot(data = data,
       mapping = aes(x = ZIE1eZiekteDt, fill = Bron)) +
  geom_bar() + 
  theme_minimal()

# epicurve of SO of hospitalised cases older than 70
ggplot(data = data %>% filter(Leeftijd >= 70),
       mapping = aes(x = ZIE1eZiekteDt, fill = Bron)) +
  geom_bar() + 
  theme_minimal()

epicurve_raw <- data %>% 
  filter(!is.na(ZIE1eZiekteDt)) %>%
  group_by(ZIE1eZiekteDt) %>%
  summarise(incidenceNL = sum(Bron != "Buitenland"),
            incidenceTOT = n()) %>%
  mutate(day = as.integer(ZIE1eZiekteDt-min(ZIE1eZiekteDt)+1),
         incidenceABROAD = incidenceTOT - incidenceNL)

epicurve_raw <- epicurve_raw %>%
  #select(day, incidenceNL, incidenceTOT) %>%
  full_join(tibble(day = 1:max(epicurve_raw$day))) %>% 
  arrange(day) %>%
  replace_na(list(incidenceNL = 0, incidenceTOT = 0, incidenceABROAD = 0)) %>%
  mutate(dates = epicurve_raw$ZIE1eZiekteDt[1] + day - 1) %>%
  select(dates, day, incidenceNL, incidenceABROAD, incidenceTOT)


#write_tsv(epicurve_raw, path = paste0("epicurve_SO_hosp_",Sys.Date(),".txt"))
#epicurve <- read_tsv("epicurve_SO_hosp_2020-03-20.txt")

# inflating reports
interval <- data$CREATIEDT - data$ZIE1eZiekteDt
interval <- interval[interval < nrow(epicurve_raw)]
reporting_delay <- ecdf(interval)
reporting_delay <- rev(reporting_delay(1:nrow(epicurve_raw)))


epicurve <- epicurve_raw %>%
  mutate(incidenceTOToriginal = incidenceTOT,
         incidenceNL = round(incidenceNL/reporting_delay),
         incidenceTOT = round(incidenceTOT/reporting_delay),
         incidenceABROAD = incidenceTOT - incidenceNL)


# use now-cast p50
epicurve <- full_join(epicurve_raw, nowcast_data, by = c("dates" = "onset.date")) %>%
  mutate(incidenceTOToriginal = incidenceTOT,
         nowcast_p50 = ifelse(is.na(nowcast_p50), incidenceTOT, nowcast_p50),
         incidenceTOT = nowcast_p50,
         incidenceABROAD = ifelse(is.na(incidenceABROAD), 0, incidenceABROAD),
         incidenceNL = incidenceTOT - incidenceABROAD,
         incidenceTOToriginal = ifelse(is.na(incidenceTOToriginal),0,incidenceTOToriginal)) 

# use now-cast p90
epicurve <- full_join(epicurve_raw, nowcast_data, by = c("dates" = "onset.date")) %>%
  mutate(incidenceTOToriginal = incidenceTOT,
         nowcast_p90 = ifelse(is.na(nowcast_p90), incidenceTOT, nowcast_p90),
         incidenceTOT = nowcast_p90,
         incidenceABROAD = ifelse(is.na(incidenceABROAD), 0, incidenceABROAD),
         incidenceNL = incidenceTOT - incidenceABROAD) 




calculate_Ru <- function(epicurve) {
  
  pdfSI <- dnorm(x = 1:7, mean = 6.8/2, sd = 1.27)
  pdfSI <- pdfSI/sum(pdfSI)
  g <- pdfSI
  
  Rt <- rep(0, nrow(epicurve))
  for(t in 2:nrow(epicurve)) {
    g2 <- g[1:min(length(g), t-1)]
    Rt[t] <- as.numeric(epicurve[t, "incidenceNL"]/sum(epicurve[t - (1:length(g2)), "incidenceTOT"]*g2))
  }
  
  Ru <- rep(0, nrow(epicurve))
  for(u in 1:(length(Rt)-1)) {
    g2 <- g[1:min(length(g), length(Rt)-u-1)]
    Ru[u] <- sum(g2*Rt[u+1:length(g2)])
  }
  
  # if(SIdelay) {
  #   SI <- rev(pnorm(q = 1:nrow(epicurve), mean = 6.8/2, sd = 1.27))
  #   Ru <- Ru/SI
  # }
  
  epicurve <- epicurve %>% mutate(Ru = Ru,
                                  Rulower = 0,
                                  Ruupper = 0)
  
  for(i in 1:nrow(epicurve)) {
    if(epicurve$incidenceTOT[i] != 0) {
      tmp <- rpois(10000, epicurve$incidenceTOT[i]*epicurve$Ru[i])/epicurve$incidenceTOT[i]
      epicurve$Rulower[i] <- quantile(tmp, probs = 0.025)
      epicurve$Ruupper[i] <- quantile(tmp, probs = 0.975)
      # tmp <- rpois(10000, epicurve$incidenceTOT[i]*epicurve$Rt[i])/epicurve$incidenceTOT[i]
      # epicurve$Rtlower[i] <- quantile(tmp, probs = 0.025)
      # epicurve$Rtupper[i] <- quantile(tmp, probs = 0.975)
    }
  }

  tmp <- which(epicurve$Ruupper == 0)
  for(i in 1:(length(tmp)-1)) {
    epicurve$Ruupper[tmp[i]] <- if(epicurve$Ruupper[tmp[i]+1] == 0) {
      (epicurve$Ruupper[tmp[i]+2] + 2*epicurve$Ruupper[tmp[i]-1])/3
    } else {
     (epicurve$Ruupper[tmp[i]+1] + epicurve$Ruupper[tmp[i]-1])/2
    }
  }
  return(epicurve)
}

plot_Ru <- function(epicurve) {
  p <- ggplot(data = epicurve,
              mapping = aes(x = dates, y = Ru, ymin = Rulower, ymax = Ruupper)) + 
    geom_ribbon(fill = "grey", col = NA) +
    geom_line() +
    geom_line(y = 1, lty = 2) +
    theme_minimal()
  return(p)
}




epicurve <- calculate_Ru(epicurve)

plot_Ru(epicurve)






# 24 maart 2020

plot_Ru <- function(epicurve) {
  par(mar = c(4,4,1,1), oma = c(0,0,0,0), las = 1, xaxs = "i", yaxs = "i")
  plot(x = epicurve$dates, y = epicurve$Ruupper, type = "l", col = "white", ylab = NA, xlab = NA)
  title(ylab = "R", line = 2)
  polygon(x = c(epicurve$dates, rev(epicurve$dates)), y = c(epicurve$Ruupper, rev(epicurve$Rulower)), border = NA, col = "grey")
  lines(x = epicurve$dates, y = epicurve$Ru)
  abline(h = 1, lty = 2)
}

p1 <- 
  ggplot(data = epicurve,
       mapping = aes(x = dates, y = incidenceTOT)) +
  geom_bar(stat = "identity") +
  geom_bar(data = epicurve,
           mapping = aes(x = dates, y = incidenceTOToriginal),
           inherit.aes = FALSE,
           stat = "identity",
           fill = 2) +
  labs(y = "aantal",
       x = "eerste ziektedag") +
  theme_minimal()

p2 <- recordPlot(plot_Ru(epicurve))

p <- plot_grid(plotlist = list(p2,p1), ncol = 1, labels = c("70+ hospitalised patients, nowcast_p50",""), label_fontface = "plain") 
ggsave(plot = p1, filename = "hospinc_24mrt2020.tiff", width = 6, height = 3, dpi = 200)



# figuur Jaap, 24 maart 2020
# met "epicurve_SO_hosp_2020-03-24.txt"
epicurve$Ruupper[1:4] <- 5

tiff(filename = "Rt_24mrt2020.tiff", width = 6, height = 3, res = 200, units = "in")
par(mar = c(4,4,1,1), oma = c(0,0,0,0), las = 1, xaxs = "i", yaxs = "i")
plot(x = epicurve$dates, y = epicurve$Ruupper, type = "l", col = "white", ylab = NA, xlab = "eerste ziektedag")
title(ylab = "R", line = 2)
polygon(x = c(epicurve$dates[1:31], rev(epicurve$dates[1:31])), y = c(epicurve$Ruupper[1:31], rev(epicurve$Rulower[1:31])), border = NA, col = "grey")
lines(x = epicurve$dates[1:31], y = epicurve$Ru[1:31])
abline(h = 1, lty = 2)
rect(xleft = as.Date("2020-02-23"), xright = as.Date("2020-02-25"), ybottom = 0.2, ytop = 0.6, col = "white")
text(x = as.Date("2020-02-24"), y = 0.4, labels = "C" )
dev.off()





# 20 maart 2020

library(EpiEstim)

incidence <- epicurve %>% 
  rename(local = incidenceNL,
         imported = incidenceABROAD) %>%
  select(dates, local, imported)
  
incidence$imported[1] <- 1
incidence$local[1] <- 0

res <- wallinga_teunis(incid = incidence, method="parametric_si",
                       config = list(t_start = seq(2, nrow(incidence)), t_end = seq(2, nrow(incidence)),
                                     mean_si = 3.6, std_si = 1.27,
                                     n_sim = 100))

plot(res)
res$R


plot(x = epicurve$dates, y = epicurve$Ruupper, type = "l", col = "white")
polygon(x = c(epicurve$dates, rev(epicurve$dates)), y = c(epicurve$Ruupper, rev(epicurve$Rulower)), border = NA, col = "grey")
lines(x = epicurve$dates, y = epicurve$Ru)
abline(h = 1, lty = 2)

plot(x = epicurve$dates, y = epicurve$Rtupper, type = "l", col = "white")
polygon(x = c(epicurve$dates, rev(epicurve$dates)), y = c(epicurve$Rtupper, rev(epicurve$Rtlower)), border = NA, col = "grey")
lines(x = epicurve$dates, y = epicurve$Rt)
abline(h = 1, lty = 2)


# Rtbootstrap <- matrix(0, nrow = 100, ncol = nrow(epicurve))
# Rubootstrap <- matrix(0, nrow = 100, ncol = nrow(epicurve))
# 
# for(i in 1:100) {
#   
#   tmp_curve <- tibble(day = base::sample(epicurve$day, size = sum(epicurve$incidenceNL), replace = TRUE, prob = epicurve$incidenceNL)) %>%
#     group_by(day) %>%
#     summarise(incidenceNL = n()) %>% 
#     full_join(tibble(day = 1:max(epicurve$day)))%>% 
#     arrange(day)%>%
#     replace_na(list(incidenceNL = 0)) %>%
#     full_join(epicurve %>% select(-incidenceNL, -incidenceTOT)) %>%
#     mutate(incidenceTOT = incidenceNL + incidenceABROAD)
#   
#   #Rt <- rep(0, nrow(tmp_curve))
#   for(t in 2:nrow(tmp_curve)) {
#     g2 <- g[1:min(length(g), t-1)]
#     Rtbootstrap[i,t] <- as.numeric(tmp_curve[t, "incidenceNL"]/sum(tmp_curve[t - (1:length(g2)), "incidenceTOT"]*g2))
#   }
#   
#   #Ru <- rep(0, nrow(tmp_curve))
#   for(u in 1:(length(Rt)-1)) {
#     g2 <- g[1:min(length(g), length(Rt)-u-1)]
#     Rubootstrap[i,u] <- sum(g2*Rt[u+1:length(g2)])
#   }
# }
# 
# Rtbootstrap[Rtbootstrap == Inf] <- NA
# Rubootstrap[Rubootstrap == Inf] <- NA
# 
# 
# boxplot(Rtbootstrap, ylim = c(0,5))
# boxplot(Rubootstrap)
# points(Ru, col = 2)
#