################################################################ ### 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 = ..." ### ### v4 = enorme verandering: tijdsafhankelijke ### overgangskansen. ### Model: ### 1. opnamebeleid in begin epidemie definieert "severe cases" (Se) => zouden zijn opgenomen aan begin van epidemie ### 2. vanuit Se is er een leeftijdsafhankelijke kans op overlijden (logistische curve) => probdeath(a) ### 3. zij die niet overlijden (C = "cured") worden allemaal opgenomen (A = "admitted") => 1 - probdeath(a) ### 4. voor alle drie tijdsafhankelijke processen die volgen, geldt ### a. tijd van een case is tijd van "severe" worden (en tijd van ziekenhuisopname indien dit gebeurt), ### ook als het gaat om processen die later plaatsvinden (ontslag, overlijden, ...) ### b. er zijn telkens vier curves als functie van tijd, voor leeftijdsgroepen 0-59, 60-69, 70-79, 80+ ### c. dit zijn logistische curves, met één parameter voor het omslagpunt (gelijk voor alle leeftijden en processen), ### één parameter voor de snelheid van de omslag (gelijk voor alle leeftijden en voor processen 2 en 3), ### en twee parameters die begin- en eindniveau bepalen (verschillend voor de vier leeftijdsgroepen en de drie processen) ### d. voor de leeftijdsgroep 0-59 wordt hierop soms een uitzondering gemaakt ### 5. tijdsafhankelijk proces 1: zij die wel overlijden (D = "dead") worden opgenomen => probzhd(a,t) ### a. het complementaire deel probdeath(a) * (1 - probzhd(a,t)) zit niet in de data ### b. voor alle leeftijdsgroepen geldt probzhd(a,t=-Inf) = 1 ### c. voor leeftijdsgroep 0-59 wordt probzhd(a,t=+Inf) = 1 verondersteld ### 6. tijdsafhankelijk proces 2: zij die overlijden en worden opgenomen gaan naar de IC => probicd(a,t) ### a. voor leeftijdsgroep 0-59 wordt probicd(a,t=+Inf) = probicd(a,t=-Inf) verondersteld (geen verandering) ### 7. tijdsafhankelijk proces 3: zij die herstellen gaan naar de IC => probicc(a,t) ### 8. er zijn nog twee tijds- en leeftijdsonafhankelijke processen: ### a. zij die overlijden, worden opgenomen en naar de IC gaan, overlijden op de IC ### (en niet op de verpleegafdeling) => probDIC ### b. zij die herstellen en naar de IC gaan, worden ontslagen vanaf de IC ### (en niet vanaf de verpleegafdeling) => probCIC ### Dit model is vervolgens omgeschreven naar overgangskansen zoals geobserveerd, en gecombineerd met kansen voor ### alle geobserveerde tijdsintervallen (delays) om de ligtijden te schatten, zodat ook de gecensorde ligtijden ### kunnen worden meegenomen ### v5: - verandering in intervallen. Afgestapt van verschillende verblijfstijden afhankelijk van uitkomst (D,C,IC/H), ### maar wel veranderende gemiddelde verblijfstijden door de tijd heen. Enige uitzondering: A2IC, want die ### is heel kort. Belangrijkste reden voor uitkomstonafhankelijk: verschillende uitkomst niet meer nodig ### voor gecensorde data (want meeste patienten zijn ontslagen), en verblijfstijd alleen gebruikt voor bezette bedden. ### Belangrijkste reden voor tijdsafhankelijk: waargenomen korte ligduur op IC. ### - extra leeftijdscategorie voor dalende mortaliteit (= dalende opnamekans), nl 50-59 ############################################# ### change points for delay distributions ### ############################################# changedates <- c(-1, 69, 147, 1000) ##################### ### 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) } analysisdate <- filedate 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) ################ ### ANALYSES ### ################ ### further clean data ### pddata <- NICEIDdataprocessed %>% filter(opnameZH %in% WhiteListedOpname) %>% mutate(ti = as.numeric(opname - as.Date("2020-02-12")), einde = if_else(inHOSP | onIC, analysisdate, einde), durA2IC = if_else(everIC, as.numeric(opnameIC - opname), NA_real_), durA2D = if_else(!everIC & HOSP2D, as.numeric(einde - opname), NA_real_), durA2C = if_else(!everIC & HOSP2C, as.numeric(einde - opname), NA_real_), durA2U = if_else(!everIC & inHOSP, as.numeric(einde - opname), NA_real_), durIC2H = if_else(everIC & (inHOSP | HOSP2D | HOSP2C), as.numeric(eindeIC - opnameIC), NA_real_), durIC2D = if_else(IC2D, as.numeric(eindeIC - opnameIC), NA_real_), durIC2C = if_else(IC2C, as.numeric(eindeIC - opnameIC), NA_real_), durIC2U = if_else(onIC, as.numeric(einde - opnameIC), NA_real_), durH2D = if_else(everIC & HOSP2D, as.numeric(einde - eindeIC), NA_real_), durH2C = if_else(everIC & HOSP2C, as.numeric(einde - eindeIC), NA_real_), durH2U = if_else(everIC & inHOSP, as.numeric(einde - eindeIC), NA_real_)) %>% mutate(ti = pmax(ti,0), age = pmin(90, age)) %>% select(ti, age, HOSP2D, HOSP2C, inHOSP, IC2D, IC2C, onIC, everIC, durA2IC, durA2D, durA2C, durA2U, durIC2H, durIC2D, durIC2C, durIC2U, durH2D, durH2C, durH2U) %>% filter(!is.na(ti) & !is.na(age)) pddata4fit <- pddata %>% mutate(age = floor(age/10), age = pmin(8, age)) ################################ ### Define all probabilities ### ################################ ### The four basic probabilities # probability to die (age) probdeath <- function(a, ad0, bd0) { exp(ad0 + bd0 * a) / (1 + exp(ad0 + bd0 * a)) } # probability to be hospitalised among those that will die probzhd <- function(a, t, pzmax, azhd, tstart) { (1 / (1 + exp(azhd * (t - tstart)))) + pzmax[a] * exp(azhd * (t - tstart)) / (1 + exp(azhd * (t - tstart))) } # probability to go to ICU among those hospitalised and about to die probicd <- function(a, t, pid0, pid1, aicd, tstart) { (pid0[a] / (1 + exp(aicd * (t - tstart)))) + pid1[a] * exp(aicd * (t - tstart)) / (1 + exp(aicd * (t - tstart))) } # probability to go to ICU among those hospitalised and about to recover probicc <- function(a, t, pic0, pic1, aicc, tstart) { (pic0[a] / (1 + exp(aicc * (t - tstart)))) + pic1[a] * exp(aicc * (t - tstart)) / (1 + exp(aicc * (t - tstart))) } ### The observation probabilities for all hospitalised # probability to die without going to ICU prob_h2d <- function(a, t, ad0, bd0, pzmax, azhd, pid0, pid1, aicd, tstart) { pdeath <- probdeath(a, ad0, bd0) pzhd <- probzhd(a, t, pzmax, azhd, tstart) picd <- probicd(a, t, pid0, pid1, aicd, tstart) return(pdeath * pzhd * (1 - picd) / (1 - pdeath + pdeath * pzhd)) } # probability to recover without going to ICU prob_h2c <- function(a, t, ad0, bd0, pzmax, azhd, pic0, pic1, aicc, tstart) { pdeath <- probdeath(a, ad0, bd0) pzhd <- probzhd(a, t, pzmax, azhd, tstart) picc <- probicc(a, t, pic0, pic1, aicc, tstart) return((1 - pdeath) * (1 - picc) / (1 - pdeath + pdeath * pzhd)) } # probability to die after going to ICU prob_i2d <- function(a, t, ad0, bd0, pzmax, azhd, pid0, pid1, aicd, tstart) { pdeath <- probdeath(a, ad0, bd0) pzhd <- probzhd(a, t, pzmax, azhd, tstart) picd <- probicc(a, t, pid0, pid1, aicd, tstart) return(pdeath * pzhd * picd / (1 - pdeath + pdeath * pzhd)) } # probability to recover after going to ICU prob_i2c <- function(a, t, ad0, bd0, pzmax, azhd, pic0, pic1, aicc, tstart) { pdeath <- probdeath(a, ad0, bd0) pzhd <- probzhd(a, t, pzmax, azhd, tstart) picc <- probicc(a, t, pic0, pic1, aicc, tstart) return((1 - pdeath) * picc / (1 - pdeath + pdeath * pzhd)) } # probability to recover prob_hi2c <- function(a, t, ad0, bd0, pzmax, azhd, tstart) { pdeath <- probdeath(a, ad0, bd0) pzhd <- probzhd(a, t, pzmax, azhd, tstart) return((1 - pdeath) / (1 - pdeath + pdeath * pzhd)) } # probability to die prob_hi2d <- function(a, t, ad0, bd0, pzmax, azhd, tstart) { pdeath <- probdeath(a, ad0, bd0) pzhd <- probzhd(a, t, pzmax, azhd, tstart) return(pdeath * pzhd / (1 - pdeath + pdeath * pzhd)) } ###################################################### ### Define all likelihoods: incremental complexity ### ###################################################### ### Likelihood 1: probability to die or recover (irrespective of ICU) # parms = {ad0, bd0, logitpzmax[5,6,7,8], logazhd, tstart, =>1-8 minloglik1 <- function(parms) { - sum(log(prob_hi2d(1 + pddata4fit$age[pddata4fit$IC2D | pddata4fit$HOSP2D], pddata4fit$ti[pddata4fit$IC2D | pddata4fit$HOSP2D], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = parms[7], tstart = parms[8]))) - sum(log(prob_hi2c(1 + pddata4fit$age[pddata4fit$IC2C | pddata4fit$HOSP2C], pddata4fit$ti[pddata4fit$IC2C | pddata4fit$HOSP2C], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = parms[7], tstart = parms[8]))) } ### Likelihood 2: L1 + probability to go to ICU among those that will die # parms = {ad0, bd0, logitpzmax[5,6,7,8], logazhd, tstart, =>1-8 # logitpid0[5,6,7,8], logitpid1[6,7,8], logaicd/logaicc, =>9-16 minloglik2 <- function(parms) { - sum(log(prob_h2d(1 + pddata4fit$age[!pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], pddata4fit$ti[!pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = exp(parms[7]), pid0 = exp(parms[1+c(8,8,8,8,8,8,9,10,11)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), pid1 = (exp(parms[1+c(8,8,8,8,8,8,12,13,14)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)]))), aicd = exp(parms[16]), tstart = parms[8]))) - sum(log(prob_i2d(1 + pddata4fit$age[pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], pddata4fit$ti[pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = exp(parms[7]), pid0 = exp(parms[1+c(8,8,8,8,8,8,9,10,11)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), pid1 = (exp(parms[1+c(8,8,8,8,8,8,12,13,14)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)]))), aicd = exp(parms[16]), tstart = parms[8]))) - sum(log(prob_hi2c(1 + pddata4fit$age[pddata4fit$IC2C | pddata4fit$HOSP2C], pddata4fit$ti[pddata4fit$IC2C | pddata4fit$HOSP2C], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = exp(parms[7]), tstart = parms[8]))) } ### Likelihood 3: L2 + probability to go to ICU among those that will recover # parms = {ad0, bd0, logitpzmax[5,6,7,8], logazhd, tstartd/tstartc, =>1-8 # logitpid0[5,6,7,8], logitpid1[6,7,8], logaicd/logaicc, =>7-16 # logitpic0[5,6,7,8], logitpic1[5,6,7,8], =>17-24 minloglik3 <- function(parms) { - sum(log(prob_h2d(1 + pddata4fit$age[!pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], pddata4fit$ti[!pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = exp(parms[7]), pid0 = exp(parms[1+c(8,8,8,8,8,8,9,10,11)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), pid1 = (exp(parms[1+c(8,8,8,8,8,8,12,13,14)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)]))), aicd = exp(parms[16]), tstart = parms[8]))) - sum(log(prob_i2d(1 + pddata4fit$age[pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], pddata4fit$ti[pddata4fit$everIC & (pddata4fit$IC2D | pddata4fit$HOSP2D)], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = exp(parms[7]), pid0 = exp(parms[1+c(8,8,8,8,8,8,9,10,11)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), pid1 = (exp(parms[1+c(8,8,8,8,8,8,12,13,14)]) / (1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)]))), aicd = exp(parms[16]), tstart = parms[8]))) - sum(log(prob_h2c(1 + pddata4fit$age[!pddata4fit$everIC & pddata4fit$HOSP2C], pddata4fit$ti[!pddata4fit$everIC & pddata4fit$HOSP2C], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = exp(parms[7]), pic0 = exp(parms[1+c(16,16,16,16,16,16,17,18,19)]) / (1 + exp(parms[1+c(16,16,16,16,16,16,17,18,19)])), pic1 = (exp(parms[1+c(20,20,20,20,20,20,21,22,23)]) / (1 + exp(parms[1+c(20,20,20,20,20,20,21,22,23)]))), aicc = exp(parms[16]), tstart = parms[8]))) - sum(log(prob_i2c(1 + pddata4fit$age[pddata4fit$everIC & (pddata4fit$HOSP2C | pddata4fit$IC2C)], pddata4fit$ti[pddata4fit$everIC & (pddata4fit$HOSP2C | pddata4fit$IC2C)], ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/ (1 + exp(parms[3:6]))), azhd = exp(parms[7]), pic0 = exp(parms[1+c(16,16,16,16,16,16,17,18,19)]) / (1 + exp(parms[1+c(16,16,16,16,16,16,17,18,19)])), pic1 = (exp(parms[1+c(20,20,20,20,20,20,21,22,23)]) / (1 + exp(parms[1+c(20,20,20,20,20,20,21,22,23)]))), aicc = exp(parms[16]), tstart = parms[8]))) } ### Likelihoods for all transitions, including delays # Admission to ICU minloglikA2IC <- function(d, A2ICr, A2ICm) { -sum(dnbinom(d, A2ICr, mu = A2ICm, log = T)) } # Admission to death minloglikA2D <- function(a, t, d, ad0, bd0, pzmax, azhd, tstart, pid0, pid1, aicd, A2Dr, A2Dm) { - sum(log(prob_h2d(1 + a, t, ad0, bd0, pzmax, azhd, pid0, pid1, aicd, tstart))) - sum(dnbinom(d, A2Dr, mu = A2Dm, log = T)) } # Admission to recovery (cure) minloglikA2C <- function(a, t, d, ad0, bd0, pzmax, azhd, tstart, pic0, pic1, aicc, A2Dr, A2Dm) { - sum(log(prob_h2c(1 + a, t, ad0, bd0, pzmax, azhd, pic0, pic1, aicc, tstart))) - sum(dnbinom(d, A2Dr, mu = A2Dm, log = T)) } # Admission but still unknown (IC not included, because that delay is very short) minloglikA2U <- function(d, A2Dr, A2Dm) { - sum(log(pnbinom(d - 1, A2Dr, mu = A2Dm, lower.tail = F))) } # ICU to regular hospital ward minloglikIC2H <- function(d, IC2Dr, IC2Dm) { -sum(dnbinom(d, IC2Dr, mu = IC2Dm, log = T)) } # ICU to death minloglikIC2D <- function(a, t, d, ad0, bd0, pzmax, azhd, tstart, pid0, pid1, aicd, IC2Dr, IC2Dm, probDIC) { - sum(log(prob_i2d(1 + a, t, ad0, bd0, pzmax, azhd, pid0, pid1, aicd, tstart))) - sum(dnbinom(d, IC2Dr, mu = IC2Dm, log = T)) - length(a) * log(probDIC) } # ICU to recovery minloglikIC2C <- function(a, t, d, ad0, bd0, pzmax, azhd, tstart, pic0, pic1, aicc, IC2Dr, IC2Dm, probCIC) { - sum(log(prob_i2c(1 + a, t, ad0, bd0, pzmax, azhd, pic0, pic1, aicc, tstart))) - sum(dnbinom(d, IC2Dr, mu = IC2Dm, log = T)) - length(a) * log(probCIC) } # ICU but still unknown minloglikIC2U <- function(d, IC2Dr, IC2Dm) { - sum(log(pnbinom(d - 1, IC2Dr, mu = IC2Dm, lower.tail = F))) } # Hospital to death minloglikH2D <- function(a, t, d, ad0, bd0, pzmax, azhd, tstart, pid0, pid1, aicd, H2Dr, H2Dm, probDIC, probreport) { - sum(log(prob_i2d(1 + a, t, ad0, bd0, pzmax, azhd, pid0, pid1, aicd, tstart))) - sum(dnbinom(d, H2Dr, mu = H2Dm, log = T)) - length(a) * (log(1 - probDIC) + log(probreport)) } # Hospital to recovery minloglikH2C <- function(a, t, d, ad0, bd0, pzmax, azhd, tstart, pic0, pic1, aicc, H2Dr, H2Dm, probCIC, probreport) { - sum(log(prob_i2c(1 + a, t, ad0, bd0, pzmax, azhd, pic0, pic1, aicc, tstart))) - sum(dnbinom(d, H2Dr, mu = H2Dm, log = T)) - length(a) * (log(1 - probCIC) + log(probreport)) } # Hospital but still unknown minloglikH2U <- function(d, H2Dr, H2Dm, probreport) { - sum(log((1 - probreport * pnbinom(d - 1, H2Dr, mu = H2Dm, lower.tail = T)))) } ### Likelihood 4: L3 + all delay distributions ### mean delays change at t = 69, t = 147 # parms = {ad0, bd0, logitpzmax[5,6,7,8], logazhd, tstartd/tstartc, =>1-8 # logitpid0[5,6,7,8], logitpid1[6,7,8], logaicd/logaicc, =>9-16 # logitpic0[5,6,7,8], logitpic1[5,6,7,8], =>17-24 # logitprobDIC, logitprobCIC, logA2ICr, logA2ICm, logA2Dr, logA2Dm[1,2,3], =>25-32 # logIC2Dr, logIC2Dm[1,2,3], logH2Dr, logH2Dm, logitprobreport =>33-39 minloglik4 <- function(parmsdelays, parmsprobs, dataseq) { parms <- c(parmsprobs, parmsdelays) minloglikA2IC(d = pddata4fit %>% filter(everIC) %>% pull(durA2IC), A2ICr = exp(parms[27]), A2ICm = exp(parms[28])) + minloglikA2D(a = pddata4fit %>% filter(!is.na(durA2D)) %>% pull(age), t = pddata4fit %>% filter(!is.na(durA2D)) %>% pull(ti), d = pddata4fit %>% filter(!is.na(durA2D)) %>% pull(durA2D), ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/(1 + exp(parms[3:6]))), azhd = exp(parms[7]), tstart = parms[8], pid0 = exp(parms[1+c(8,8,8,8,8,8,9,10,11)])/(1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), pid1 = exp(parms[1+c(8,8,8,8,8,8,12,13,14)])/(1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)])), aicd = exp(parms[16]), A2Dr = exp(parms[29]), A2Dm = exp(parms[30:32])[ pddata4fit %>% filter(!is.na(durA2D)) %>% mutate(per = sapply(ti, function(x) sum(x > dataseq))) %>% pull(per) ] ) + minloglikA2C(a = pddata4fit %>% filter(!is.na(durA2C)) %>% pull(age), t = pddata4fit %>% filter(!is.na(durA2C)) %>% pull(ti), d = pddata4fit %>% filter(!is.na(durA2C)) %>% pull(durA2C), ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/(1 + exp(parms[3:6]))), azhd = exp(parms[7]), tstart = parms[8], pic0 = exp(parms[1+c(16,16,16,16,16,16,17,18,19)])/(1 + exp(parms[1+c(16,16,16,16,16,16,17,18,19)])), pic1 = exp(parms[1+c(20,20,20,20,20,20,21,22,23)])/(1 + exp(parms[1+c(20,20,20,20,20,20,21,22,23)])), aicc = exp(parms[16]), A2Dr = exp(parms[29]), A2Dm = exp(parms[30:32])[ pddata4fit %>% filter(!is.na(durA2C)) %>% mutate(per = sapply(ti, function(x) sum(x > dataseq))) %>% pull(per) ] ) + minloglikA2U(d = pddata4fit %>% filter(!is.na(durA2U)) %>% pull(durA2U), A2Dr = exp(parms[29]), A2Dm = exp(parms[30:32])[ pddata4fit %>% filter(!is.na(durA2U)) %>% mutate(per = sapply(ti, function(x) sum(x > dataseq))) %>% pull(per) ] ) + minloglikIC2H(d = pddata4fit %>% filter(!is.na(durIC2H)) %>% pull(durIC2H), IC2Dr = exp(parms[33]), IC2Dm = exp(parms[34:36])[ pddata4fit %>% filter(!is.na(durIC2H)) %>% mutate(per = sapply(ti, function(x) sum(x > dataseq))) %>% pull(per) ] ) + minloglikIC2D(a = pddata4fit %>% filter(!is.na(durIC2D)) %>% pull(age), t = pddata4fit %>% filter(!is.na(durIC2D)) %>% pull(ti), d = pddata4fit %>% filter(!is.na(durIC2D)) %>% pull(durIC2D), ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/(1 + exp(parms[3:6]))), azhd = exp(parms[7]), tstart = parms[8], pid0 = exp(parms[1+c(8,8,8,8,8,8,9,10,11)])/(1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), pid1 = exp(parms[1+c(8,8,8,8,8,8,12,13,14)])/(1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)])), aicd = exp(parms[16]), IC2Dr = exp(parms[33]), IC2Dm = exp(parms[34:36])[ pddata4fit %>% filter(!is.na(durIC2D)) %>% mutate(per = sapply(ti, function(x) sum(x > dataseq))) %>% pull(per) ], probDIC = exp(parms[25])/(1 + exp(parms[25]))) + minloglikIC2C(a = pddata4fit %>% filter(!is.na(durIC2C)) %>% pull(age), t = pddata4fit %>% filter(!is.na(durIC2C)) %>% pull(ti), d = pddata4fit %>% filter(!is.na(durIC2C)) %>% pull(durIC2C), ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/(1 + exp(parms[3:6]))), azhd = exp(parms[7]), tstart = parms[8], pic0 = exp(parms[1+c(16,16,16,16,16,16,17,18,19)])/(1 + exp(parms[1+c(16,16,16,16,16,16,17,18,19)])), pic1 = exp(parms[1+c(20,20,20,20,20,20,21,22,23)])/(1 + exp(parms[1+c(20,20,20,20,20,20,21,22,23)])), aicc = exp(parms[16]), IC2Dr = exp(parms[33]), IC2Dm = exp(parms[34:36])[ pddata4fit %>% filter(!is.na(durIC2C)) %>% mutate(per = sapply(ti, function(x) sum(x > dataseq))) %>% pull(per) ], probCIC = exp(parms[26])/(1 + exp(parms[26]))) + minloglikIC2U(d = pddata4fit %>% filter(!is.na(durIC2U)) %>% pull(durIC2U), IC2Dr = exp(parms[33]), IC2Dm = exp(parms[34:36])[ pddata4fit %>% filter(!is.na(durIC2U)) %>% mutate(per = sapply(ti, function(x) sum(x > dataseq))) %>% pull(per) ] ) + minloglikH2D(a = pddata4fit %>% filter(!is.na(durH2D)) %>% pull(age), t = pddata4fit %>% filter(!is.na(durH2D)) %>% pull(ti), d = pddata4fit %>% filter(!is.na(durH2D)) %>% pull(durH2D), ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/(1 + exp(parms[3:6]))), azhd = exp(parms[7]), tstart = parms[8], pid0 = exp(parms[1+c(8,8,8,8,8,8,9,10,11)])/(1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), pid1 = exp(parms[1+c(8,8,8,8,8,8,12,13,14)])/(1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)])), aicd = exp(parms[16]), H2Dr = exp(parms[37]), H2Dm = exp(parms[38]), probDIC = exp(parms[25])/(1 + exp(parms[25])), probreport = exp(parms[39])/(1 + exp(parms[39]))) + minloglikH2C(a = pddata4fit %>% filter(!is.na(durH2C)) %>% pull(age), t = pddata4fit %>% filter(!is.na(durH2C)) %>% pull(ti), d = pddata4fit %>% filter(!is.na(durH2C)) %>% pull(durH2C), ad0 = parms[1], bd0 = parms[2], pzmax = c(rep(1, 5), exp(parms[3:6])/(1 + exp(parms[3:6]))), azhd = exp(parms[7]), tstart = parms[8], pic0 = exp(parms[1+c(16,16,16,16,16,16,17,18,19)])/(1 + exp(parms[1+c(16,16,16,16,16,16,17,18,19)])), pic1 = exp(parms[1+c(20,20,20,20,20,20,21,22,23)])/(1 + exp(parms[1+c(20,20,20,20,20,20,21,22,23)])), aicc = exp(parms[16]), H2Dr = exp(parms[37]), H2Dm = exp(parms[38]), probCIC = exp(parms[26])/(1 + exp(parms[26])), probreport = exp(parms[39])/(1 + exp(parms[39]))) + minloglikH2U(d = pddata4fit %>% filter(!is.na(durH2U)) %>% pull(durH2U), H2Dr = exp(parms[37]), H2Dm = exp(parms[38]), probreport = exp(parms[39])/(1 + exp(parms[39]))) } ### step by step estimation: initial values to improve convergence; warnings may still occur parms1 <- optim(c(-8,1,.5,0.5,0.5,-1,0.5,40), minloglik1, method = "BFGS")$par parms2 <- optim(c(parms1, c(0.5,2,0.5,-2,0,-1.5,-4,-2.5)), minloglik2, method = "BFGS")$par parms3 <- optim(c(parms2, c(-0.5,-0.5,-0.5,-4,-2,-1.5,-3,-5)), minloglik3, method = "BFGS")$par # optim(c(2,2), # function(x) minloglik4(c(parms3, x, rep(2,13))), method = "BFGS") parms4 <- optim(c(2,-5,-1,.5,.5,2,2,2,0.5,3,3,2.5,1,2,4.5), minloglik4, parmsprobs = parms3, dataseq = changedates, method = "BFGS")$par NICEparameters <- c(parms3, parms4) ### check with latest results that required adaptation (data upto 10 sept 2020) if(any(abs(NICEparameters - c(-7.1,0.8,-0.7,-0.3,0.1,-0.8,-2.1,42.6, 0.4,1.1,0.0,-2.1,0.2,-0.9,-3.2,-1.8, -0.8,-0.6,-0.9,-3.8,-1.7,-1.3,-1.9,-4.0, 2.3,-4.7,-1.0,.7,0.1,2.2,2.3,2.1,0.3,3,2.7,2.6,0.6,2.3,4.9)) > 0.5)) { print("Please check NICE results carefully") } ######################################################## ### Functions to obtain all probabilities and delays ### ######################################################## ### function to calculate the probabilities get_allprobs <- function(parms, tmax) { a <- matrix(rep(1:9, each = tmax + 1), ncol = 9) t <- matrix(rep(0:tmax, 9), ncol = 9) pD <- probdeath(a, parms[1], parms[2]) pAinD <- probzhd(a, t, c(rep(1, 5), exp(parms[3:6])/(1 + exp(parms[3:6]))), exp(parms[7]), parms[8]) pICinAD <- probicd(a, t, exp(parms[1+c(8,8,8,8,8,8,9,10,11)])/(1 + exp(parms[1+c(8,8,8,8,8,8,9,10,11)])), exp(parms[1+c(8,8,8,8,8,8,12,13,14)])/(1 + exp(parms[1+c(8,8,8,8,8,8,12,13,14)])), exp(parms[16]), parms[8]) pDICinICAD <- exp(parms[25])/(1 + exp(parms[25])) pICinC <- probicc(a, t, exp(parms[1+c(16,16,16,16,16,16,17,18,19)])/(1 + exp(parms[1+c(16,16,16,16,16,16,17,18,19)])), exp(parms[1+c(20,20,20,20,20,20,21,22,23)])/(1 + exp(parms[1+c(20,20,20,20,20,20,21,22,23)])), exp(parms[16]), parms[8]) pCICinICC <- exp(parms[26])/(1 + exp(parms[26])) probSe2A <- pD * pAinD + 1 - pD probSe2D <- pD * (1 - pAinD) probA2IC <- (pD * pAinD * pICinAD + (1 - pD) * pICinC) / probSe2A probA2D <- pD * pAinD * (1 - pICinAD) / probSe2A probA2C <- (1 - pD) * (1 - pICinC) / probSe2A probIC2H <- (pD * pAinD * pICinAD * (1 - pDICinICAD) + (1 - pD) * pICinC * (1 - pCICinICC)) / probSe2A / probA2IC probIC2D <- pD * pAinD * pICinAD * pDICinICAD / probSe2A / probA2IC probIC2C <- (1 - pD) * pICinC * pCICinICC / probSe2A / probA2IC probH2D <- pD * pAinD * pICinAD * (1 - pDICinICAD) / probSe2A / probA2IC / probIC2H probH2C <- (1 - pD) * pICinC * (1 - pCICinICC) / probSe2A / probA2IC / probIC2H probSe2A2IC <- pD * pAinD * pICinAD + (1 - pD) * pICinC probSe2A2D <- pD * pAinD * (1 - pICinAD) probSe2A2C <- (1 - pD) * (1 - pICinC) probSe2A2IC2H <- (pD * pAinD * pICinAD * (1 - pDICinICAD) + (1 - pD) * pICinC * (1 - pCICinICC)) probSe2A2IC2D <- pD * pAinD * pICinAD * pDICinICAD probSe2A2IC2C <- (1 - pD) * pICinC * pCICinICC probSe2A2IC2H2D <- pD * pAinD * pICinAD * (1 - pDICinICAD) probSe2A2IC2H2C <- (1 - pD) * pICinC * (1 - pCICinICC) return(list( probSe2A = probSe2A, probSe2D = probSe2D, probA2IC = probA2IC, probA2D = probA2D, probA2C = probA2C, probIC2H = probIC2H, probIC2D = probIC2D, probIC2C = probIC2C, probH2D = probH2D, probH2C = probH2C, probSe2A2IC = probSe2A2IC, probSe2A2D = probSe2A2D, probSe2A2C = probSe2A2C, probSe2A2IC2H = probSe2A2IC2H, probSe2A2IC2D = probSe2A2IC2D, probSe2A2IC2C = probSe2A2IC2C, probSe2A2IC2H2D = probSe2A2IC2H2D, probSe2A2IC2H2C = probSe2A2IC2H2C )) } ### function to calculate the delays get_alldelays <- function(parms, maxdelay = 100) { delayA2IC <- dnbinom(0:maxdelay, exp(parms[27]), mu = exp(parms[28])) delayA2D <- sapply(exp(parms[30:32]), function(x) dnbinom(0:maxdelay, exp(parms[29]), mu = x)) delayIC2D <- sapply(exp(parms[34:36]), function(x) dnbinom(0:maxdelay, exp(parms[33]), mu = x)) delayH2D <- dnbinom(0:maxdelay, exp(parms[37]), mu = exp(parms[38])) delayA2IC2D <- t(sapply(1:(1 + maxdelay), function(x) colSums(delayA2IC[1:x] * delayIC2D[x:1,,drop = F]))) delayA2IC2H2D <- t(sapply(1:(1 + maxdelay), function(x) colSums(delayA2IC2D[1:x,,drop=F] * delayH2D[x:1]))) return(list( delayA2IC = delayA2IC, delayA2D = delayA2D, delayIC2D = delayIC2D, delayH2D = delayH2D, delayA2IC2D = delayA2IC2D, delayA2IC2H2D = delayA2IC2H2D, changedates = changedates )) } ######################## ### 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, maxligduur = Inf) { 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, eindeIC) %>% 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")), rehospdag = as.numeric(eindeIC - as.Date("2020-02-12")), rehospdag = if_else(is.na(rehospdag), einddag, rehospdag), einddag = pmin(einddag, rehospdag + maxligduur)) %>% select(hospdag, einddag) %>% filter() 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))) } } } ################################################## ### Create lists with probabilities and delays ### ################################################## NICEprobabilities <- get_allprobs(NICEparameters, as.numeric(analysisdate - as.Date("2020-02-12"))) NICEdelays <- get_alldelays(NICEparameters, as.numeric(analysisdate - as.Date("2020-02-12"))) ################################ ### Remove temporary objects ### ################################ rm(parms1, parms2, parms3, probdeath, probzhd, probicc, probicd, prob_h2c, prob_h2d, prob_i2c, prob_i2d, prob_hi2c, prob_hi2d, minloglikA2C, minloglikA2D, minloglikA2IC, minloglikA2U, minloglikIC2C, minloglikIC2D, minloglikIC2H, minloglikIC2U, minloglikH2C, minloglikH2D, minloglikH2U, minloglik1, minloglik2, minloglik3, minloglik4, changedates)