not_ready/benchmark_test_mellem_drg.R

In davidbaniadam/farver: Some Useful Function

# se den anden fil for en "samlet benchmark" for de forskellige procestider.
#### Forskel i middel indenfor forskellige operationsgrupper med konfidensintervalller.


# library(tidyverse)
# library(lubridate)
# library(readxl)
# library(infer)
# library(stringr)
# library(rlang)

# load data ---------------------------------------------------------------
# Det er denne datasæt koden arbejder med. Ser kun på Urologi (én speciale)
load("L:/LovbeskyttetMapper/BogA - Analyse/DBAN/Benchmark(unden kobling)/RH3_og_HEH3_OP_data.RData")

# To do liste:
# Rens data for na?
# evt. med median
# boksplot med outliers
# ggplot så title kan komme med.






# Først skal vi finde faktisk forskel (i samlettid).
# Vi sætter de to datasæt ovenpå hinadnen og beregner "samlet_tid".
# Her beregner vi gns. for hver drg grupperet på Enhed (dvs. hospital).
# vi forskelle i gns for hver Enhed (hospital) for forskellige drg.
# Og så tager vi differencen mellem forskellige enheder.

# Laver bootstrap fordeling for "forskel i middel" over alle drg grupper.
# Beregner konfidensintervaller på 3 forskellige måder. I plotten er Bais-corrrected KI med.
# Skal laves en funktion så reps i boot og i formlen for Bais-corrected KI kan være en variabel.

forskel <- bind_rows(OP_aktivitet_H3, OP_aktivitet_R3) %>%
  mutate(foertid = as.numeric(difftime(Knivtid_Proc_Start, Patient_Paa_Opstue, units="mins")),
         knivtid = as.numeric(difftime(Sidste_Sutur_Proc_Slut, Knivtid_Proc_Start, units="mins")),
         eftertid = as.numeric(difftime(Patient_Forlader_Stuen, Sidste_Sutur_Proc_Slut, units="mins")),
         samlet_tid = as.numeric(difftime(Patient_Forlader_Stuen, Patient_Paa_Opstue, units="mins")))  %>%
  group_by(drg) %>%
  nest()

rm(OP_aktivitet_H3, OP_aktivitet_R3)



bootloop <- function(dataset, procestid, method=1, reps = 1000, alpha = 0.05) {
  procestid <- enquo(procestid)

  diff_mean <- dataset %>%
    mutate(diff_means  = map(data, function(.x){.x %>%
        group_by(Enhed) %>%
        summarise(mean(!!procestid, na.rm=TRUE)) %>%
        pull() %>%
        diff() })) %>%
    select(-data) # Pass på at diff passer med nedenstående i "calculate".


  bootstrap <- dataset %>%
    mutate(distribution =map(data, function(.x){ .x %>%
        specify(as.formula(paste0(quo_name(procestid), "~ Enhed")) ) %>%
        generate(reps = reps, type = "bootstrap") %>%
        calculate(stat = "diff in means", order = c( "RH Urologi", "HeH Urologisk Kirurgi H"))} )) %>% # NB: pas på at det er samme diff som diff_mean
    inner_join(diff_mean, by="drg")

  if (method==1){
    bootstrap2 <- bootstrap %>% mutate(Bias_Corrected_KI=map2(distribution, diff_means, function(.x, .y){ .x %>%
        summarise( l =quantile(.x$stat,pnorm(2*qnorm(sum(.x$stat >= .y)/reps) + qnorm(alpha/2))),
                   u= quantile(.x$stat,pnorm(2*qnorm(sum(.x$stat >= .y)/reps) + qnorm(1-alpha/2)))    )})) }# Her bliver der lavet Bais-corrected KI

  else if (method==2){
    bootstrap2 <- bootstrap %>% mutate(Percentile_KI = map(distribution, function(.x){.x %>%
        summarize(l = quantile(stat, alpha/2),
                  u = quantile(stat, 1 - alpha/2))})) }
  else{
    bootstrap2 <- bootstrap %>% mutate(SD_KI =map2(distribution, diff_means, function(.x,.y){.x %>%
        conf_int(level = (1 - alpha), type="se", point_estimate = .y)}))
  }

  bootstrap3 <- bootstrap2 %>% mutate(plots= map2(distribution, Bias_Corrected_KI, function(.x, .y){ .x  %>%
      visualize(bins=20) +  shade_confidence_interval(.y)
  }))
  return(bootstrap3)
}


aa <-bootloop(forskel, samlet_tid, method=1)
aa$plots
procestider <- list("foertid", "knivtid", "eftertid", "samlet_tid")


a <- map(syms(procestider), bootloop , dataset=forskel, reps=1000)

warnings()

a[[4]]$Bias_Corrected_KI

# warnings pga. visualize.
# Brug ggplot i sted for visualize så titel kan komme med.
# ggplot(diff_mean_ci, aes(stat)) + geom_histogram(fill="blue", bins = 20) +
# ggtitle("Konfidensinterval for forskel i middel") + geom_vline(xintercept=diff_mean)