################################################################
### Functions to estimate probabilities and delays from NICE ###
################################################################

###   v1 = correcte behandeling van IC-hosp-transitie:  ###
###           vanuit IC zijn er drie paden:             ###
###           dood, ontslag, ZH. Het ZH-pad werd fout   ###
###           gecodeerd, nl met de status ná het ZH     ###
###           (dood of ontslag). Nu correct.            ###

###   v2 = onbetrouwbare ziekenhuizen worden      ###
###   eruit gehaald voor kansen en delays         ###
###   vanuit het ziekenhuis                       ###

###   v3 = typo in maken NICEIDdataprocessed:         ###
###   haakje voor na.rm = T ipv erna in "einde = ..." ###


#####################
### Read the data ###
#####################

### read NICE file with IDs
files <- list.files("data", full.names = TRUE)
files <- files[grepl(files, pattern = "NICEID")]
filedates <- as.Date(substr(files, 13, 20), "%d%m%Y")
if(exists("analysisdate")) {
  file.name <- files[filedates == analysisdate]
  filedate <- analysisdate
} else {
  file.name <- files[filedates %>% which.max]
  filedate <- max(filedates)
}
NICEIDdata <- read_csv(file.name, 
                       col_types = cols(
                         pid = col_double(),
                         seq = col_double(),
                         name = col_character(),
                         age = col_character(),
                         adm_date_icu = col_character(),
                         dis_date_icu = col_character(),
                         discharged_to = col_character(),
                         is_from_other_ic = col_character(),
                         naar_zkh = col_character(),
                         is_ic = col_double(),
                         died_in_hospital = col_character(),
                         covid19status = col_character()
                       )) 
NICEIDdataprocessed <- NICEIDdata %>%
  mutate(age = if_else(age != "null", as.numeric(age), NA_real_),
         adm_date_icu = as.Date(adm_date_icu, "%d/%m/%Y"),
         dis_date_icu = as.Date(dis_date_icu, "%d/%m/%Y"),
         discharged_to = case_when(
           discharged_to == "null" ~ "censored",
           discharged_to %in% as.character(c(1, 3, 4, 6)) ~ "HOSP",
           discharged_to %in% as.character(c(8)) ~ "home",
           discharged_to == "7" ~ "death",
           discharged_to == "9" ~ "unknown",
           discharged_to %in% as.character(c(2, 5, 10)) ~ "otherICU",
           is_ic == 1 ~ "censored",
           TRUE ~ "hospitalrecord"
         )) %>%
  arrange(pid, seq) %>%
  group_by(pid) %>% 
  # only patients that were tested positive
  filter(any(covid19status %in% c("lab", "ct"))) %>%
  mutate(
    # sometimes age is different in different hospitals (typo)
    age = max(age),
    # first hospital that is not mentioned in 'referred to'
    opnameZH = setdiff(name, naar_zkh)[1],
    # solution in case all hospitals were referred to (patient went back and forth)
    opnameZH = if_else(any(is.na(opnameZH)), name[seq == min(seq)], opnameZH[1]),
    laatsteZH = tail(name[adm_date_icu == max(adm_date_icu)], 1)
  ) %>% 
  mutate(
    # first admitted at IC/HOSP
    admitted_at = if_else(any(is_ic[adm_date_icu == min(c(adm_date_icu), na.rm = T)] == 1), 
                          "IC", "HOSP", "UNK"),
    everIC = any(is_ic == 1),
    movementIC2HOSP = if_else(
      (any(is_ic == 1) & any(is_ic == 0) & 
         min(c(filedate + 100, adm_date_icu[is_ic == 1]), na.rm = T) < 
         max(c(filedate - 1000, adm_date_icu[is_ic == 0]), na.rm = T)) |
        any(discharged_to == "HOSP"),
      TRUE, FALSE, FALSE
    ),
    movementHOSP2IC = if_else(
      any(is_ic == 1) & any(is_ic == 0) & 
        min(c(filedate + 100, adm_date_icu[is_ic == 0]), na.rm = T) < 
        max(c(filedate - 1000, adm_date_icu[is_ic == 1]), na.rm = T),
      TRUE, FALSE, FALSE
    ),
    IC2D = any(discharged_to == "death"),
    IC2C = any(discharged_to == "home"),
    IC2G = any(discharged_to == "Germany"),
    HOSP2D = any(died_in_hospital == 1) & !IC2D,
    HOSP2C = all(!is.na(dis_date_icu)) & !IC2D & !IC2C& !HOSP2D,
    onIC = !IC2D & !IC2C & !IC2G  & !HOSP2D & !HOSP2C & 
      any(is_ic[adm_date_icu == max(adm_date_icu)] == 1 &
            discharged_to[adm_date_icu == max(adm_date_icu)] == "censored"),
    inHOSP = !IC2D & !IC2C & !HOSP2D & !HOSP2C & !onIC,
    opname = min(adm_date_icu),
    opnameIC = if_else(any(is_ic == 1), 
                       min(c(filedate + 100, adm_date_icu[is_ic == 1])), 
                       as.Date(NA)),
    einde = if_else(IC2D | IC2C  | IC2G | HOSP2D | HOSP2C,
                    max(c(filedate - 1000, dis_date_icu), na.rm = T), 
                    as.Date(NA)),
    eindeIC = if_else(any(is_ic == 1) & !onIC, 
                      max(c(filedate - 1000, dis_date_icu[is_ic == 1])), 
                      as.Date(NA)),
    eindeIC = if_else(everIC & !onIC & is.na(eindeIC) & 
                        max(c(filedate - 1000, adm_date_icu)) > max(c(filedate - 1000, adm_date_icu[is_ic == 1])),
                      max(adm_date_icu, na.rm = T),
                      eindeIC,
                      eindeIC)
  ) %>%
  ungroup() %>%
  select(pid, age, opnameZH, laatsteZH, admitted_at, everIC, movementIC2HOSP,
         movementHOSP2IC, IC2D, IC2C, IC2G, HOSP2D, HOSP2C, onIC, inHOSP,
         opname, opnameIC, einde, eindeIC) %>%
  distinct()


WhiteListedOpname <- NICEIDdataprocessed %>%
  group_by(opnameZH) %>%
  summarise(pctIC = sum(everIC)/n()) %>%
  filter(pctIC < 0.6) %>% pull(opnameZH)
WhiteListedOntslag <- NICEIDdataprocessed %>%
  group_by(laatsteZH) %>%
  filter(any(!everIC)) %>%
  summarise(pctinhosp = sum(inHOSP[!everIC])/sum(!everIC)) %>% 
  filter(pctinhosp < 0.25) %>% pull(laatsteZH)



#######################
### hospital to ICU ###
#######################

get_H2IC_delays <- function() {
  NICEIDdataprocessed %>%
    filter(everIC & opnameZH %in% WhiteListedOpname) %>%
    mutate(interval = as.numeric(opnameIC - opname)) %>%
    select(interval) %>% 
    summarise(meanH2IC = mean(interval),
              varH2IC = var(interval),
              pH2IC = 1 - meanH2IC/varH2IC,
              rH2IC = meanH2IC * (1-pH2IC)/pH2IC,
              totaalH2ICd = n()) %>%
    select(meanH2IC, pH2IC, rH2IC, totaalH2ICd) %>%
    as.list()
}


get_H2IC_probs <- function() {
  NICEIDdataprocessed %>%
    filter(opnameZH %in% WhiteListedOpname) %>%
    filter(!is.na(age)) %>%
    mutate(
      ageclass = if_else(is.na(age), -1, floor(age/10)),
      ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
    group_by(ageclass) %>%
    summarise(totaalH2ICp = n(),
              H2IC = sum(everIC),
              probH2IC = H2IC/totaalH2ICp) %>%
    select(totaalH2ICp, H2IC, probH2IC) %>%
    as.list()
}


##############################
### Hospital to death/cure ###
##############################

HOSPintervals <- NICEIDdataprocessed %>%
  filter(opname > as.Date("2020-02-12")) %>%
  filter(!everIC) %>%
  filter(laatsteZH %in% WhiteListedOntslag) %>%
  mutate(ageclass = floor(age / 10),
         ageclass = pmin(ageclass, 8),
         death = IC2D | HOSP2D,
         cured = IC2C | HOSP2C,
         censored = !death & !cured,
         ligtijd = if_else(death | cured,
                           as.numeric(einde - opname),
                           as.numeric(filedate - opname))) %>%
  filter(!is.na(ligtijd))


minlogLikH2DC <- function(pars) {
  pdeath <- exp(pars[5] + HOSPintervals$ageclass * pars[6])/
    (1 + exp(pars[5] + HOSPintervals$ageclass * pars[6]))
  pdeath[is.na(pdeath)] <- mean(pdeath, na.rm = T)

  shD <- exp(pars[1])
  muD <- exp(pars[2])
  shC <- exp(pars[3])
  muC <- exp(pars[4])
  
  - sum(log(pdeath[HOSPintervals$death] * 
              dnbinom(HOSPintervals$ligtijd[HOSPintervals$death], 
                      shD, mu = muD))) -
    sum(log((1 - pdeath[HOSPintervals$cured]) * 
              dnbinom(HOSPintervals$ligtijd[HOSPintervals$cured], 
                      shC, mu = muC))) -
    sum(log(pdeath[HOSPintervals$censored] * 
              pnbinom(HOSPintervals$ligtijd[HOSPintervals$censored] - 1, 
                      shD, mu = muD, lower.tail = FALSE) +
              (1 - pdeath[HOSPintervals$censored]) * 
              pnbinom(HOSPintervals$ligtijd[HOSPintervals$censored] - 1, 
                      shC, mu = muC, lower.tail = FALSE)))
}

get_H2DC_probsdelays <- function() {
  res <- optim(c(1,1,1,1,0,0), minlogLikH2DC)
  return(
    list(
      probH2D = exp(res$par[5] + (0:8) * res$par[6]) / (1 + exp(res$par[5] + (0:8) * res$par[6])),
      meanH2D = exp(res$par[2]),
      rH2D = exp(res$par[1]),
      pH2D = exp(res$par[2])/(exp(res$par[2]) + exp(res$par[1])),
      meanH2C = exp(res$par[4]),
      rH2C = exp(res$par[3]),
      pH2C = exp(res$par[4])/(exp(res$par[4]) + exp(res$par[3]))
    )
  )
}


##################################
### ICU to hospital/death/cure ###
##################################

ICintervals <- NICEIDdataprocessed %>%
  filter(opname > as.Date("2020-02-12")) %>%
  filter(everIC) %>%
  mutate(ageclass = floor(age / 10),
         death = IC2D | HOSP2D,
         cured = IC2C | HOSP2C,
         tohosp = inHOSP | HOSP2D | HOSP2C,
         censored = onIC | IC2G,
         ligtijd = if_else(death | cured | tohosp | IC2G,
                           as.numeric(eindeIC - opnameIC),
                           as.numeric(filedate - opnameIC))) %>%
  filter(!is.na(ligtijd))

minlogLikIC2DHC <- function(pars) {
  pdeath <- exp(pars[7] + ICintervals$ageclass * pars[8])/
    (1 + exp(pars[7] + ICintervals$ageclass * pars[8]))
  pdeath[is.na(pdeath)] <- mean(pdeath, na.rm = T)
  
  phosp <- exp(pars[9])/ (1 + exp(pars[9]))
  phosp <- phosp * (1 - pdeath)
  

  shD <- exp(pars[1])
  muD <- exp(pars[2])
  shH <- exp(pars[3])
  muH <- exp(pars[4])
  shC <- exp(pars[5])
  muC <- exp(pars[6])
  
  - sum(log(pdeath[ICintervals$IC2D] * 
              dnbinom(ICintervals$ligtijd[ICintervals$IC2D], 
                      shD, mu = muD))) -
    sum(log(phosp[ICintervals$tohosp] *
              dnbinom(ICintervals$ligtijd[ICintervals$tohosp],
                      shH, mu = muH))) -
    sum(log((1 - pdeath[ICintervals$IC2C] - phosp[ICintervals$IC2C]) * 
              dnbinom(ICintervals$ligtijd[ICintervals$IC2C], 
                      shC, mu = muC))) -
    sum(log(pdeath[ICintervals$censored] * 
              pnbinom(ICintervals$ligtijd[ICintervals$censored] - 1, 
                      shD, mu = muD, lower.tail = FALSE) +
              phosp[ICintervals$censored] * 
              pnbinom(ICintervals$ligtijd[ICintervals$censored] - 1, 
                      shH, mu = muH, lower.tail = FALSE) +
              (1 - pdeath[ICintervals$censored] - phosp[ICintervals$censored]) * 
              pnbinom(ICintervals$ligtijd[ICintervals$censored] - 1, 
                      shC, mu = muC, lower.tail = FALSE)))
}

get_IC2DHC_probsdelays <- function() {
  res <- optim(c(.5,2.5,.5,3,0,2,-5,0,1), minlogLikIC2DHC)
  return(
    list(
      probIC2D = exp(res$par[7] + (0:8) * res$par[8]) / (1 + exp(res$par[7] + (0:8) * res$par[8])),
      probIC2H = exp(res$par[9]) / (1 + exp(res$par[9])),
      meanIC2D = exp(res$par[2]),
      rIC2D = exp(res$par[1]),
      pIC2D = exp(res$par[2])/(exp(res$par[2]) + exp(res$par[1])),
      meanIC2H = exp(res$par[4]),
      rIC2H = exp(res$par[3]),
      pIC2H = exp(res$par[4])/(exp(res$par[4]) + exp(res$par[3])),
      meanIC2C = exp(res$par[6]),
      rIC2C = exp(res$par[5]),
      pIC2C = exp(res$par[6])/(exp(res$par[6]) + exp(res$par[5]))
    )
  )
}


##################################
### 2nd hospital to death/cure ###
##################################
HOSP2intervals <- NICEIDdataprocessed %>%
  filter(opname > as.Date("2020-02-12")) %>%
  filter(everIC & !is.na(eindeIC) & (is.na(einde) | einde > eindeIC)) %>%
  filter(laatsteZH %in% WhiteListedOntslag) %>%
  mutate(ageclass = floor(age / 10),
         death = IC2D | HOSP2D,
         cured = IC2C | HOSP2C,
         censored = !death & !cured,
         ligtijd = if_else(death | cured,
                           as.numeric(einde - eindeIC),
                           as.numeric(filedate - eindeIC))) %>%
  filter(!is.na(ligtijd))

# same distribution for both outcomes until more data available
minlogLikH22DC <- function(pars) {
  pdeath <- exp(pars[5] + HOSP2intervals$ageclass * pars[6])/
    (1 + exp(pars[5] + HOSP2intervals$ageclass * pars[6]))
  pdeath[is.na(pdeath)] <- mean(pdeath, na.rm = T)
  
  shD <- exp(pars[1])
  muD <- exp(pars[2])
  shC <- exp(pars[1])
  muC <- exp(pars[2])
  
  - sum(log(pdeath[HOSP2intervals$death] * 
              dnbinom(HOSP2intervals$ligtijd[HOSP2intervals$death], 
                      shD, mu = muD))) -
    sum(log((1 - pdeath[HOSP2intervals$cured]) * 
              dnbinom(HOSP2intervals$ligtijd[HOSP2intervals$cured], 
                      shC, mu = muC))) -
    sum(log(pdeath[HOSP2intervals$censored] * 
              pnbinom(HOSP2intervals$ligtijd[HOSP2intervals$censored] - 1, 
                      shD, mu = muD, lower.tail = FALSE) +
              (1 - pdeath[HOSP2intervals$censored]) * 
              pnbinom(HOSP2intervals$ligtijd[HOSP2intervals$censored] - 1, 
                      shC, mu = muC, lower.tail = FALSE)))
}

get_H22DC_probsdelays <- function() {
  res <- optim(c(.6,2.6,.6,2.6,0,0), minlogLikH22DC)
  return(
    list(
      probH22D = exp(res$par[5] + (0:8) * res$par[6]) / (1 + exp(res$par[5] + (0:8) * res$par[6])),
      meanH22D = exp(res$par[2]),
      rH22D = exp(res$par[1]),
      pH22D = exp(res$par[2])/(exp(res$par[2]) + exp(res$par[1])),
      meanH22C = exp(res$par[2]),
      rH22C = exp(res$par[1]),
      pH22C = exp(res$par[2])/(exp(res$par[2]) + exp(res$par[1]))
    )
  )
}


########################
### incidence curves ###
########################

inc_nice_hosp <- function(ages = FALSE, ROAZ = "all") {
  if(is.numeric(ROAZ)) {
    ROAZ <- sort(unique(ROAZdata$Regio_ROAZ))[-9][ROAZ]
  }
  if(ROAZ == "all") {
    if(ages) {
      HOSPdata <- NICEIDdataprocessed %>%
        select(age, opname) %>%
        filter(opname <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(HOSPdag = as.numeric(opname - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(HOSPdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(HOSPdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        select(opname) %>%
        filter(opname <= filedate) %>%
        mutate(hospdag = as.numeric(opname - as.Date("2020-02-12"))) %>%
        pull(hospdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  } else {
    if(ages) {
      HOSPdata <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(age, opname) %>%
        filter(opname <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(HOSPdag = as.numeric(opname - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(HOSPdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(HOSPdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(opname) %>%
        filter(opname <= filedate) %>%
        mutate(hospdag = as.numeric(opname - as.Date("2020-02-12"))) %>%
        pull(hospdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  }
}


inc_nice_ic <- function(ages = FALSE, ROAZ = "all") {
  if(is.numeric(ROAZ)) {
    ROAZ <- sort(unique(ROAZdata$Regio_ROAZ))[-9][ROAZ]
  }
  if(ROAZ == "all") {
    if(ages) {
      ICdata <- NICEIDdataprocessed %>%
        select(age, opnameIC) %>%
        filter(opnameIC <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(ICdag = as.numeric(opnameIC - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(ICdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(ICdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        select(opnameIC) %>%
        filter(opnameIC <= filedate) %>%
        mutate(icdag = as.numeric(opnameIC - as.Date("2020-02-12"))) %>%
        pull(icdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  } else {
    if(ages) {
      ICdata <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(age, opnameIC) %>%
        filter(opnameIC <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(ICdag = as.numeric(opnameIC - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(ICdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(ICdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(opnameIC) %>%
        filter(opnameIC <= filedate) %>%
        mutate(icdag = as.numeric(opnameIC - as.Date("2020-02-12"))) %>%
        pull(icdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  }
}

dis_nice_hosp <- function(ages = FALSE, ROAZ = "all") {
  if(is.numeric(ROAZ)) {
    ROAZ <- sort(unique(ROAZdata$Regio_ROAZ))[-9][ROAZ]
  }
  if(ROAZ == "all") {
    if(ages) {
      HOSPdata <- NICEIDdataprocessed %>%
        select(age, einde) %>%
        filter(einde <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(HOSPdag = as.numeric(einde - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(HOSPdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(HOSPdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        select(einde) %>%
        filter(!is.na(einde) & einde <= filedate) %>%
        mutate(hospdag = as.numeric(einde - as.Date("2020-02-12"))) %>%
        pull(hospdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  } else {
    if(ages) {
      HOSPdata <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(age, einde) %>%
        filter(einde <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(HOSPdag = as.numeric(einde - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(HOSPdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(HOSPdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(einde) %>%
        filter(!is.na(einde) & einde <= filedate) %>%
        mutate(hospdag = as.numeric(einde - as.Date("2020-02-12"))) %>%
        pull(hospdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  }
}

dis_nice_ic <- function(ages = FALSE, ROAZ = "all") {
  if(is.numeric(ROAZ)) {
    ROAZ <- sort(unique(ROAZdata$Regio_ROAZ))[-9][ROAZ]
  }
  if(ROAZ == "all") {
    if(ages) {
      ICdata <- NICEIDdataprocessed %>%
        select(age, eindeIC) %>%
        filter(eindeIC <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(ICdag = as.numeric(eindeIC - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(ICdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(ICdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        select(eindeIC) %>%
        filter(!is.na(eindeIC) & eindeIC <= filedate) %>%
        mutate(icdag = as.numeric(eindeIC - as.Date("2020-02-12"))) %>%
        pull(icdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  } else {
    if(ages) {
      ICdata <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(age, eindeIC) %>%
        filter(eindeIC <= filedate) %>%
        filter(!is.na(age)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass)) %>%
        group_by(ageclass) %>%
        mutate(ICdag = as.numeric(eindeIC - as.Date("2020-02-12"))) 
      toreturn <- lapply(0:8, function(x) tabulate(ICdata %>% 
                                                     filter(ageclass == x) %>% 
                                                     pull(ICdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
      
    } else {
      NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(eindeIC) %>%
        filter(!is.na(eindeIC) & eindeIC <= filedate) %>%
        mutate(icdag = as.numeric(eindeIC - as.Date("2020-02-12"))) %>%
        pull(icdag) %>% tabulate(nbins = as.numeric(filedate - as.Date("2020-02-12")))
    }
  }
}

#########################
### prevalence curves ###
#########################

prev_nice_hosp <- function(ages = FALSE, ROAZ = "all", corrected = FALSE, incdelay = NA, disdelay = NA) {
  if(is.numeric(ROAZ)) {
    ROAZ <- sort(unique(ROAZdata$Regio_ROAZ))[-9][ROAZ]
  }
  if(ROAZ == "all") {
    if(ages) {
      mutations <- NICEIDdataprocessed %>%
        select(opname, einde, age) %>%
        filter(opname <= filedate & (is.na(einde) | einde <= filedate)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass),
          einde = if_else(is.na(einde), filedate + 1, einde)) %>%
        mutate(hospdag = as.numeric(opname - as.Date("2020-02-12")),
               einddag = as.numeric(einde - as.Date("2020-02-12")))
      opnames <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                    filter(ageclass == x) %>% 
                                                    pull(hospdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      ontslagen <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                      filter(ageclass == x) %>% 
                                                      pull(einddag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      if(corrected) {
        opnames <- lapply(opnames, function(x) x/rev(incdelay[1:length(x)]))
        ontslagen <- lapply(ontslagen, function(x) x/rev(disdelay[1:length(x)]))
      }
      toreturn <- lapply(1:9, function(x) round(cumsum(opnames[[x]]) - cumsum(ontslagen[[x]])))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
    } else {
      mutations <- NICEIDdataprocessed %>%
        select(opname, einde) %>%
        filter(opname <= filedate & (is.na(einde) | einde <= filedate)) %>%
        mutate(hospdag = as.numeric(opname - as.Date("2020-02-12")),
               einde = if_else(is.na(einde), filedate + 1, einde),
               einddag = as.numeric(einde - as.Date("2020-02-12"))) %>%
        select(hospdag, einddag)
      opnames <- tabulate(mutations$hospdag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      ontslagen <- tabulate(mutations$einddag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      if(corrected) {
        opnames <- opnames / rev(incdelay[1:length(opnames)])
        ontslagen <- ontslagen / rev(disdelay[1:length(ontslagen)])
      }
      return(round(cumsum(opnames) - cumsum(ontslagen)))
    }
  } else {
    if(ages) {
      mutations <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(opname, einde, age) %>%
        filter(opname <= filedate & (is.na(einde) | einde <= filedate)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass),
          einde = if_else(is.na(einde), filedate + 1, einde)) %>%
        mutate(hospdag = as.numeric(opname - as.Date("2020-02-12")),
               einddag = as.numeric(einde - as.Date("2020-02-12")))
      opnames <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                    filter(ageclass == x) %>% 
                                                    pull(hospdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      ontslagen <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                      filter(ageclass == x) %>% 
                                                      pull(einddag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      if(corrected) {
        opnames <- lapply(opnames, function(x) x/rev(incdelay[1:length(x)]))
        ontslagen <- lapply(ontslagen, function(x) x/rev(disdelay[1:length(x)]))
      }
      toreturn <- lapply(1:9, function(x) round(cumsum(opnames[[x]]) - cumsum(ontslagen[[x]])))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
    } else {
      mutations <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(opname, einde) %>%
        filter(opname <= filedate & (is.na(einde) | einde <= filedate)) %>%
        mutate(hospdag = as.numeric(opname - as.Date("2020-02-12")),
               einde = if_else(is.na(einde), filedate + 1, einde),
               einddag = as.numeric(einde - as.Date("2020-02-12"))) %>%
        select(hospdag, einddag)
      opnames <- tabulate(mutations$hospdag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      ontslagen <- tabulate(mutations$einddag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      if(corrected) {
        opnames <- opnames / rev(incdelay[1:length(opnames)])
        ontslagen <- ontslagen / rev(disdelay[1:length(ontslagen)])
      }
      return(round(cumsum(opnames) - cumsum(ontslagen)))
    }
  }
}

prev_nice_ic <- function(ages = FALSE, ROAZ = "all", corrected = FALSE, incdelay = NA, disdelay = NA) {
  if(is.numeric(ROAZ)) {
    ROAZ <- sort(unique(ROAZdata$Regio_ROAZ))[-9][ROAZ]
  }
  if(ROAZ == "all") {
    if(ages) {
      mutations <- NICEIDdataprocessed %>%
        select(opnameIC, eindeIC, age) %>%
        filter(opnameIC <= filedate & (is.na(eindeIC) | eindeIC <= filedate)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass),
          eindeIC = if_else(is.na(eindeIC), filedate + 1, eindeIC)) %>%
        mutate(ICdag = as.numeric(opnameIC - as.Date("2020-02-12")),
               einddag = as.numeric(eindeIC - as.Date("2020-02-12")))
      opnames <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                    filter(ageclass == x) %>% 
                                                    pull(ICdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      ontslagen <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                      filter(ageclass == x) %>% 
                                                      pull(einddag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      if(corrected) {
        opnames <- lapply(opnames, function(x) x/rev(incdelay[1:length(x)]))
        ontslagen <- lapply(ontslagen, function(x) x/rev(disdelay[1:length(x)]))
      }
      toreturn <- lapply(1:9, function(x) round(cumsum(opnames[[x]]) - cumsum(ontslagen[[x]])))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
    } else {
      mutations <- NICEIDdataprocessed %>%
        select(opnameIC, eindeIC) %>%
        filter(opnameIC <= filedate & (is.na(eindeIC) | eindeIC <= filedate)) %>%
        mutate(icdag = as.numeric(opnameIC - as.Date("2020-02-12")),
               eindeIC = if_else(is.na(eindeIC), filedate + 1, eindeIC),
               einddag = as.numeric(eindeIC - as.Date("2020-02-12"))) %>%
        select(icdag, einddag)
      opnames <- tabulate(mutations$icdag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      ontslagen <- tabulate(mutations$einddag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      if(corrected) {
        opnames <- opnames / rev(incdelay[1:length(opnames)])
        ontslagen <- ontslagen / rev(disdelay[1:length(ontslagen)])
      }
      return(round(cumsum(opnames) - cumsum(ontslagen)))
    }
  } else {
    if(ages) {
      mutations <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(opnameIC, eindeIC, age) %>%
        filter(opnameIC <= filedate & (is.na(eindeIC) | eindeIC <= filedate)) %>%
        mutate(
          ageclass = floor(age / 10),
          ageclass = if_else(ageclass > 8, 8, ageclass),
          eindeIC = if_else(is.na(eindeIC), filedate + 1, eindeIC)) %>%
        mutate(ICdag = as.numeric(opnameIC - as.Date("2020-02-12")),
               einddag = as.numeric(eindeIC - as.Date("2020-02-12")))
      opnames <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                    filter(ageclass == x) %>% 
                                                    pull(ICdag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      ontslagen <- lapply(0:8, function(x) tabulate(mutations %>% 
                                                      filter(ageclass == x) %>% 
                                                      pull(einddag), nbins = as.numeric(filedate - as.Date("2020-02-12"))))
      if(corrected) {
        opnames <- lapply(opnames, function(x) x/rev(incdelay[1:length(x)]))
        ontslagen <- lapply(ontslagen, function(x) x/rev(disdelay[1:length(x)]))
      }
      toreturn <- lapply(1:9, function(x) round(cumsum(opnames[[x]]) - cumsum(ontslagen[[x]])))
      names(toreturn) <- c("[0,10)", "[10,20)", "[20,30)", "[30,40)",
                           "[40,50)", "[50,60)", "[60,70)", "[70,80)", "[80,Inf]")
      return(toreturn)
    } else {
      mutations <- NICEIDdataprocessed %>%
        filter(ZH_ROAZ == ROAZ) %>%
        select(opnameIC, eindeIC) %>%
        filter(opnameIC <= filedate & (is.na(eindeIC) | eindeIC <= filedate)) %>%
        mutate(icdag = as.numeric(opnameIC - as.Date("2020-02-12")),
               einde = if_else(is.na(eindeIC), filedate + 1, eindeIC),
               einddag = as.numeric(eindeIC - as.Date("2020-02-12"))) %>%
        select(icdag, einddag)
      opnames <- tabulate(mutations$icdag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      ontslagen <- tabulate(mutations$einddag, nbins = as.numeric(filedate - as.Date("2020-02-12")))
      if(corrected) {
        opnames <- opnames / rev(incdelay[1:length(opnames)])
        ontslagen <- ontslagen / rev(disdelay[1:length(ontslagen)])
      }
      return(round(cumsum(opnames) - cumsum(ontslagen)))
    }
  }
}