farver: Some Useful Function

# Benchmark overordnet.
# Definitioner og forkortelser brugt:
# PH: primær hospital, BM: benchmark hospital.

# Kort Beskrivelse af kode:
# i- Denne kode tager to datasæt med operationstider som input og laver en benchmark analyse.
# ii- Ser kun på én speciale.
# iii- drg grupperne er tilfældige tal i  koden  og det skal justeres senere.
# iv- ser  på forskel i procent i procestiderne mellem de to hospitaler fx "samlet tid".
# v- Konfidensintervallerne bliver lavet via bootstrap-metoden vha pakken "infer".





# Ting som skal tilføjes senere:
# kontrollere for komobilitet evt. via brug af "covariate?"; brug køn og alder til at estimere...
# Evt. loopes så der laves benchmark på flere hospitaler.
# Udskiftningstid skal beregnes for sig selv hvis vi skal lave bootstrap.
library(tidyverse)
library(lubridate)
# library(readxl)
# library(infer)


# Først beregner vi den faktiske forskel i procent mellem PH og BM

BM <- OP_aktivitet_R3 %>%
  group_by(drg) %>%
  summarise(BM_mean_samlet_tid = mean(as.numeric(difftime(Patient_Forlader_Stuen, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            BM_mean_foertid =    mean(as.numeric(difftime(Knivtid_Proc_Start, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            BM_mean_knivtid =    mean(as.numeric(difftime(Sidste_Sutur_Proc_Slut, Knivtid_Proc_Start, units="mins")), na.rm=TRUE),
            BM_mean_eftertid =   mean(as.numeric(difftime(Patient_Forlader_Stuen, Sidste_Sutur_Proc_Slut, units="mins")), na.rm=TRUE),
            BM_antal  =     n())



PH <- OP_aktivitet_H3 %>%
  group_by(drg) %>%
  summarise(PH_mean_samlet_tid = mean(as.numeric(difftime(Patient_Forlader_Stuen, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            PH_mean_foertid =    mean(as.numeric(difftime(Knivtid_Proc_Start, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            PH_mean_knivtid =    mean(as.numeric(difftime(Sidste_Sutur_Proc_Slut, Knivtid_Proc_Start, units="mins")), na.rm=TRUE),
            PH_mean_eftertid =   mean(as.numeric(difftime(Patient_Forlader_Stuen, Sidste_Sutur_Proc_Slut, units="mins")), na.rm=TRUE),
            PH_antal  =     n())


# Samler BH og PH og udregner faktisk forskel i procent som er givet ved PH/BM -1
# Faktiske forskelle i procent.
# Vi bergner total antal minuter brugt på operation i PH og BM (i BM holder vi bruger vi "antal" fra PH)
samlet <- inner_join(BM, PH, by = c( "drg" = "drg")) %>%
  mutate(total_PH_mean_samlet_tid    = PH_antal*PH_mean_samlet_tid,
         total_BM_mean_samlet_tid_if = PH_antal*BM_mean_samlet_tid,

         total_PH_mean_foertid      = PH_antal*PH_mean_foertid,
         total_BM_mean_foertid_if   = PH_antal*BM_mean_foertid,

         total_PH_mean_knivtid      = PH_antal*PH_mean_knivtid,
         total_BM_mean_knivtid_if   = PH_antal*BM_mean_knivtid,

         total_PH_mean_eftertid     = PH_antal*PH_mean_eftertid,
         total_BM_mean_eftertid_if  = PH_antal*BM_mean_eftertid) %>%

  summarise(samlet_tid_forskel_i_pct =sum(total_PH_mean_samlet_tid) / sum(total_BM_mean_samlet_tid_if)-1,

            foertid_forskel_i_pct     =sum(total_PH_mean_foertid)/sum(total_BM_mean_foertid_if)-1,

            knivtid_forskel_i_pct   =sum(total_PH_mean_knivtid)/sum(total_BM_mean_knivtid_if)-1,

            eftertid_forskel_i_pct   =sum(total_PH_mean_eftertid)/sum(total_BM_mean_eftertid_if)-1)

# Her laver vi konfidensintervaller for forskellene i procent

BM_CI <- OP_aktivitet_R3 %>%
  generate(reps=1000, type="bootstrap") %>%
  group_by(drg, replicate) %>%
  summarise(BM_mean_samlet_tid = mean(as.numeric(difftime(Patient_Forlader_Stuen, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            BM_mean_foertid =    mean(as.numeric(difftime(Knivtid_Proc_Start, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            BM_mean_knivtid =    mean(as.numeric(difftime(Sidste_Sutur_Proc_Slut, Knivtid_Proc_Start, units="mins")), na.rm=TRUE),
            BM_mean_eftertid =   mean(as.numeric(difftime(Patient_Forlader_Stuen, Sidste_Sutur_Proc_Slut, units="mins")), na.rm=TRUE),
            BM_antal  =     n())



PH_CI <- OP_aktivitet_H3 %>%
  generate(reps=1000, type="bootstrap") %>%
  group_by(drg, replicate) %>%
  summarise(PH_mean_samlet_tid = mean(as.numeric(difftime(Patient_Forlader_Stuen, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            PH_mean_foertid =    mean(as.numeric(difftime(Knivtid_Proc_Start, Patient_Paa_Opstue, units="mins")), na.rm=TRUE),
            PH_mean_knivtid =    mean(as.numeric(difftime(Sidste_Sutur_Proc_Slut, Knivtid_Proc_Start, units="mins")), na.rm=TRUE),
            PH_mean_eftertid =   mean(as.numeric(difftime(Patient_Forlader_Stuen, Sidste_Sutur_Proc_Slut, units="mins")), na.rm=TRUE),
            PH_antal  =     n())



# Med denne kode kan vi finde K.I. for forskellige procestider givet i "forskel i procent".

samlet_CI_join <- inner_join(BM_CI, PH_CI, by = c("replicate" = "replicate", "drg" = "drg")) %>%
  mutate(total_PH_mean_samlet_tid    = PH_antal*PH_mean_samlet_tid,
         total_BM_mean_samlet_tid_if = PH_antal*BM_mean_samlet_tid,

         total_PH_mean_foertid      = PH_antal*PH_mean_foertid,
         total_BM_mean_foertid_if   = PH_antal*BM_mean_foertid,

         total_PH_mean_knivtid      = PH_antal*PH_mean_knivtid,
         total_BM_mean_knivtid_if   = PH_antal*BM_mean_knivtid,

         total_PH_mean_eftertid     = PH_antal*PH_mean_eftertid,
         total_BM_mean_eftertid_if  = PH_antal*BM_mean_eftertid) %>%

  group_by(replicate) %>%
  summarise(samlet_tid_forskel_i_pct =sum(total_PH_mean_samlet_tid) / sum(total_BM_mean_samlet_tid_if)-1,

            foertid_forskel_i_pct     =sum(total_PH_mean_foertid)/sum(total_BM_mean_foertid_if)-1,

            knivtid_forskel_i_pct   =sum(total_PH_mean_knivtid)/sum(total_BM_mean_knivtid_if)-1,

            eftertid_forskel_i_pct   =sum(total_PH_mean_eftertid)/sum(total_BM_mean_eftertid_if)-1)


# Metode 1: Konfidensinterval metoden
# 95% konfidensinterval via " Percentile method"
#l for lower og u for upper
samlet_CI_join %>%
  summarise(l_samlet_tid = quantile(samlet_tid_forskel_i_pct, 0.025),
            u_samlet_tid = quantile(samlet_tid_forskel_i_pct, 0.975),

            l_foertid = quantile(foertid_forskel_i_pct, 0.025),
            u_foertid = quantile(foertid_forskel_i_pct, 0.975),

            l_knivtid = quantile(knivtid_forskel_i_pct, 0.025),
            u_knivtid = quantile(knivtid_forskel_i_pct, 0.975),

            l_eftertid = quantile(eftertid_forskel_i_pct, 0.025),
            u_eftertid = quantile(eftertid_forskel_i_pct, 0.975))


# 95% konfidensinterval via sd metoden.
#  Man kan også bruge  qt(0.975, df = ?) isted for 1.96 for små stikprøver.
samlet_CI_join %>%
  summarise(boot_se = sd(samlet_tid_forskel_i_pct)) %>%
  summarise(l = pull(samlet[1,1]) - 1.96 * boot_se,
            u = pull(samlet[1,1]) + 1.96 * boot_se)



# Histogram over bootstrap fordelingen for "samlet tid" og vertikal line som viser selve estimationen.
ggplot(samlet_CI_join, aes(samlet_tid_forskel_i_pct)) + geom_histogram(fill="blue", bins = 20) +  geom_vline(xintercept=pull(samlet[1,1])) +
  ggtitle("Fordeling over forskel i pct")