# pop <- haven::read_sas(here::here("E:/rawdata/709203/Population/pop_talos.sas7bdat"))

# sst <- list.files(here::here("E:/rawdata/709203/Eksterne data"), pattern = "*.sas7bdat", full.names = TRUE) |>
#   purrr::map(haven::read_sas)

source(here::here("R/glmnet-reg.R"))

#' Read all sas files in folder to list
#'
#' @param path folder path
#'
#' @return list
sas2list <- function(path) {
  ls <- list.files(here::here(path), pattern = "*.sas7bdat", full.names = TRUE) |>
    purrr::map(haven::read_sas)
  names(ls) <- list.files(here::here(path), pattern = "*.sas7bdat") |>
    gsub(".sas7bdat", "", x = _) |>
    toupper()
  ls
}


#' Flatten multilevel list
#'
#' @param paths character vector of folder paths
#'
#' @return flattened list
flatmultiread <- function(paths) {
  paths |>
    purrr::map(sas2list) |>
    purrr::list_flatten()
}

# ls <- targets::tar_read(reg_data)

# Vectors are kept for compatibility. Calling functions can be done from within other functions. So much easier!

date_cutter <- function() as.Date("2023-01-01")
date.cut <- date_cutter() # The earliest date will define the overall date cut

censor_cutter <- function() 12
censor.cut <- censor_cutter() # years of maximum follow up, due to small numbers

vasc.diags <- c("I21", "I61", "I63", "I64","G45", "K28")

#' Extract deaths from Dødsårsagsregiseret
#'
#' @param ls 
#' @param max.date 
#' @param diags.vasc 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(reg_data) |> get_deaths()
get_deaths <- function(ls, max.date = date.cut, diags.vasc = vasc.diags) {
  
  
  # vasc.death.tilg <- mapply(ls$DAR_T_DODSAARSAG_2[, "C_DODTILGRUNDL_ACME"],
  #   FUN = function(i) { # mapply inside apply call to handle rowvise matching in in matrix
  #     i_3 <- substr(i, 1, 3) # Substr only the 3 first characters to match by group
  #     i_3 %in% diags.vasc # Rowvise matching
  #   }
  # )
  
  vasc.death.tilg <- ls$DAR_T_DODSAARSAG_2[["C_DODTILGRUNDL_ACME"]] |> substr(1, 3) %in% diags.vasc
  

  vasc.death.any <- apply(mapply(ls$DAR_T_DODSAARSAG_2[, c("C_DODTILGRUNDL_ACME", "C_DOD_1A", "C_DOD_1B", "C_DOD_1C", "C_DOD_1D")],
    FUN = function(i) { # mapply inside apply call to handle rowvise matching in in matrix
      i_3 <- substr(i, 1, 3) # Substr only the 3 first characters to match by group
      i_3 %in% diags.vasc # Rowvise matching
    }
  ), 1, any) # Simplify to TRUE if any
  #
  vasc.death.other <- apply(mapply(ls$DAR_T_DODSAARSAG_2[, c("C_DOD_1A", "C_DOD_1B", "C_DOD_1C", "C_DOD_1D")],
    FUN = function(i) { # mapply inside apply call to handle rowvise matching in in matrix
      i_3 <- substr(i, 1, 3) # Substr only the 3 first characters to match by group
      i_3 %in% diags.vasc # Rowvise matching
    }
  ), 1, any) # Simplify to TRUE if any

  diag.either <- xor(vasc.death.other, vasc.death.tilg)
  diag.both <- vasc.death.other & vasc.death.tilg

  vasc.death.diags <-
    apply(
      mapply(
        ls$DAR_T_DODSAARSAG_2[, c(
          "C_DODTILGRUNDL_ACME",
          "C_DOD_1A",
          "C_DOD_1B",
          "C_DOD_1C",
          "C_DOD_1D"
        )],
        FUN = function(i) {
          # mapply inside apply call to handle rowvise matching in in matrix
          substr(i, 1, 3) # Substr only the 3 first characters to match by group
        }
      ),
      1,
      paste,
      collapse = ","
    )

  deaths.vasc <-
    ls$DAR_T_DODSAARSAG_2 |>
    dplyr::select(K_CPR, D_STATDATO) |>
    dplyr::filter(vasc.death.tilg)

  df.death.all <- ls$CPR3_T_PERSON |>
    dplyr::filter(C_STATUS == 90) |> # People migrating are filtered (n ~ 1)
    dplyr::select(c(
      "V_PNR",
      "D_STATUS_HEN_START"
    )) |>
    dplyr::left_join(deaths.vasc, by = c("V_PNR" = "K_CPR")) |>
    dplyr::mutate(vasc_death = !is.na(D_STATDATO)) |>
    dplyr::transmute(
      PNR = V_PNR,
      death_date = D_STATUS_HEN_START,
      vasc_death = vasc_death
    )

  df.death.all |> dplyr::filter(death_date < max.date)
}


#' Title
#'
#' @param ls 
#' @param max.date 
#' @param diags.vasc 
#'
#' @return
#' @export
#'
#' @examples
#' ls <- targets::tar_read(reg_list)
get_events <- function(ls, max.date = date.cut, diags.vasc = vasc.diags) {
  ident.vars <- toupper(c("_recnum", "_cpr"))

  df.vasc.events.lpr <- ls$LPR_T_DIAG |>
    dplyr::mutate(dia.f = substr(C_DIAG, 2, 4)) |> # Subsets only 2:4 chars, to get main group
    dplyr::filter(
      dia.f %in% diags.vasc
      #     &  # Filters to only include pre-defined diagnoses
      #       C_DIAGTYPE=="A"
    ) |> # Filters to only include if main diagnosis
    dplyr::select(
      ends_with(ident.vars),
      "C_DIAG",
      "C_DIAGTYPE"
    ) |>
    dplyr::left_join(
      ls$LPR_T_ADM |> dplyr::select(
        tidyselect::ends_with(ident.vars),
        "D_INDDTO",
        "C_INDM",
        "D_UDDTO",
        "C_UDM",
        "C_SGH",
        "C_AFD",
        "C_ADIAG"
      ),
      by = c("V_RECNUM" = "K_RECNUM")
    )

  ## LPR-F - LPR 3

  ident.vars.lpr3 <- toupper(c("cpr", "_kontakt", "DW_EK_FORLOEB"))

  df.vasc.events.lpr3 <- ls$LPR_F_DIAGNOSER |>
    dplyr::mutate(dia.f = substr(DIAGNOSEKODE, 2, 4)) |> # Subsets only 2:4 chars, to get main group
    dplyr::filter(
      dia.f %in% diags.vasc
      #       &  # Filters to only include pre-defined diagnoses
      #         DIAGNOSETYPE=="A"
    ) |> # Filters to only include if main diagnosis
    dplyr::select(
      ends_with(ident.vars.lpr3),
      "DIAGNOSEKODE",
      "DIAGNOSETYPE"
    ) |>
    dplyr::left_join(
      ls$LPR_F_KONTAKTER |> dplyr::select(
        ends_with(ident.vars.lpr3),
        "DATO_START",
        "DATO_SLUT",
        "PRIORITET"
      ),
      by = c("DW_EK_KONTAKT")
    ) |>
    dplyr::mutate(PRIORITET = as.character((PRIORITET == "ATA1") + 1)) # If ATA1 then 1, if not (ATA3) then 2, cowboy coding

  df.vasc.events <- dplyr::full_join(df.vasc.events.lpr, df.vasc.events.lpr3, by = c(
    "V_CPR" = "CPR",
    "C_DIAG" = "DIAGNOSEKODE",
    "C_DIAGTYPE" = "DIAGNOSETYPE",
    "D_INDDTO" = "DATO_START",
    "D_UDDTO" = "DATO_SLUT",
    "C_INDM" = "PRIORITET"
  )) |>
    dplyr::mutate(date.event = D_INDDTO)

  df.vasc.events |> dplyr::filter(date.event < max.date)
}


# deaths <- targets::tar_read(df_deaths)
# events <- targets::tar_read(df_events)
# clinical <- targets::tar_read(pop_df)


#' Filter only truly considered events
#'
#' @param data
#'
#' @return tibble
define_events <- function(data) {
  data |> dplyr::filter(
    C_DIAGTYPE == "A", # Primary diagnosis
    C_INDM == "1", # Acutely admitted
    difftime(date.event, rdate, units = "days") > 5 # More than five (5) days after randomisation/primary stroke
  )
}

#' The big merger and filter of events
#'
#' @param ls list of events, deaths and clinical
#'
#' @return tibble
#' 
#' @examples
#' ls <- list(events = targets::tar_read(df_events), deaths = targets::tar_read(df_deaths), clinical = targets::tar_read(pop_df))
#' ls |> all_events()
all_events <- function(ls) {
  # ls <- list(events = targets::tar_read(df_events), deaths = targets::tar_read(df_deaths), clinical = targets::tar_read(pop_df))
  df.events <- dplyr::full_join(
    purrr::pluck(ls, "events"),
    purrr::pluck(ls, "deaths") |>
      dplyr::mutate(
        # These are just added to ease later filtering
        C_DIAGTYPE = "A",
        C_INDM = "1"
      ), # Ads diagtype=A, C_INDM=1 for easier sorting later
    by = c(
      "V_CPR" = "PNR",
      "date.event" = "death_date",
      "C_DIAGTYPE",
      "C_INDM"
    )
  ) |>
    dplyr::full_join(dplyr::select(purrr::pluck(ls, "clinical"), c("PNR", "rdate", "enddate")),
      by = c("V_CPR" = "PNR")
    ) |>
    dplyr::arrange(date.event) |> # Sort by event date
    dplyr::mutate(event.type = dplyr::if_else(
      is.na(C_DIAG),
      dplyr::if_else(vasc_death, "death.vasc", "death.other"),
      substr(C_DIAG, 1, 4)
    ))

  df.events |>
    dplyr::group_split(V_CPR) |> # Splits by CPR
    purrr::map(define_events) |> # Custom function to specify criteria for events
    purrr::discard(\(x) nrow(x) == 0) |>   # Discard empty elements
    purrr::map(dplyr::transmute, # Saving only relevant variables
               CPR = V_CPR,
               date.event,
               event.type)
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' ls <- list(events = targets::tar_read(df_events), deaths = targets::tar_read(df_deaths), clinical = targets::tar_read(pop_df))
#' ls |> merge_events()
merge_events <- function(data){
  data|>
  all_events() |>
    purrr::modify(\(x) x[1, ]) |> # Select first event
    purrr::list_rbind()
}


#' Count number of prescriptions of given ATC group for each CPR
#'
#' @param atc.code atc group
#' @param data dataset from LMS
#'
#' @return tibble
count_treat <- function(atc.code, data) {
  data |>
    dplyr::filter(grepl(atc.code, ATC)) |>
    dplyr::count(CPR) |>
    dplyr::filter(n > 1)
}


#' Get LMS data
#'
#' @param ls list of datasets
#' @param max.date filter date for max inclusion
#' @param atc.tbl tibble of atc codes
#'
#' @return tibble
get_lms <- function(ls, max.date, atc.tbl) {
  ls$LMS_EPIKUR |>
    dplyr::filter(grepl(atc.tbl[1], ATC)) |>
    dplyr::left_join(ls$LMS_LAEGEMIDDELOPLYSNINGER) |>
    dplyr::filter(as.Date(ACTDATE) < max.date)
}

#' Get count of treated patients from LMS data
#'
#' @param ls list of datasets
#' @param max.date filter date for max inclusion
#' @param atc.tbl tibble of atc codes
#'
#' @return tibble
get_treated <- function(ls,
                        max.date = date.cut,
                        atc.tbl = c(
                          atc.antidep = "N06A",
                          atc.ssri = "N06AB"
                        )) {
  df <- atc.tbl |>
    purrr::map(count_treat, get_lms(ls, max.date, atc.tbl)) |>
    purrr::reduce(dplyr::full_join, by = "CPR")
  colnames(df) <- c("CPR", paste0("n.", names(atc.tbl)))
  df
}

#' BMI calc, drops
#'
#' @param w weight in kg
#' @param h height in cm
#' @param data data set
#'
#' @return tibble
bmi_calc <- function(data, drop = TRUE) {
  # After inspection, both h+w are missing if any is missing
  out <- data |> dplyr::mutate(reg_bmi = suppressWarnings(as.numeric(dplyr::if_else(reg_vaegt == "NA", reg_vaegt_anslaaet, reg_vaegt)) / ((as.numeric(reg_hojde) / 100)^2)))

  if (drop) {
    out <- out |>
      dplyr::select(-dplyr::all_of(c("reg_vaegt", "reg_vaegt_anslaaet", "reg_hojde")))
  }
  out
}

is_equal <- function(data, test) {
  data == test
}

#' Load clinical population data
#'
#' @return tibble
#' @examples
#' get_clinical() |> colnames()
#' 
get_clinical <- function() {
  sas2list("E:/rawdata/709203/Population")[[2]] |> 
    correct_na() |>
    dplyr::mutate(
      reg_smoker = dplyr::case_match(
        reg_rygning, "1" ~ TRUE,
        c("2", "3", "4") ~ FALSE,
        "9" ~ NA
      ),
      # Living alone defined as not together with somebody
      reg_alone = dplyr::case_match(
        reg_civil, "1" ~ FALSE,
        c("2", "3") ~ TRUE,
        "9" ~ NA
      ),
      reg_more_alc = dplyr::case_match(
        reg_alkohol, "1" ~ FALSE,
        "2" ~ TRUE,
        "9" ~ NA
      ),
      reg_female = sex == "Kvinde",
      dplyr::across(
        .cols = c(
          "reg_hyperten",
          "reg_diabetes",
          "reg_atriefli",
          "reg_perifer_arteriel",
          "reg_tidl_tci",
          "reg_ami"
        ),
        ~ dplyr::case_match(
          .x, "1" ~ TRUE,
          "2" ~ FALSE,
          "9" ~ NA
        )
      ),
      dplyr::across(
        .cols = c(
          "reg_trombolyse",
          "reg_trombektomi"
        ),
        ~ dplyr::case_match(
          .x, "1" ~ TRUE,
          c("3","4") ~ FALSE,
          "9" ~ NA
        )
      ),
      reg_any_perf = reg_trombolyse | reg_trombektomi
    ) |>
    bmi_calc()
}

# get_clinical() |> pragmatic_imputation() |> skimr::skim()

# get_clinical <- function() {
#   sas2list("E:/rawdata/709203/Population")[[2]] |>
#     dplyr::mutate(
#       reg_smoker = reg_rygning == 1,
#       reg_cohabiting = reg_civil == 1,
#       reg_more_alc = reg_alkohol == 2,
#       reg_female = sex == "Kvinde",
#       dplyr::across(.cols = c("reg_hyperten", "reg_diabetes", "reg_atriefli", "reg_perifer_arteriel", "reg_tidl_tci", "reg_ami", "reg_trombolyse", "reg_trombektomi"), ~ .x == 1),
#       reg_any_perf = reg_trombolyse | reg_trombektomi
#     ) |>
#     bmi_calc()
# }

pragmatic_imputation <- function(data, vec=c("reg_trombolyse", "reg_trombektomi","reg_any_perf")) {
  # Assumes, if not TRUE, then FALSE (gets rid of NAs)
  data |> dplyr::mutate(dplyr::across(.cols = tidyselect::any_of(vec), ~dplyr::if_else(.x,TRUE,FALSE,missing = FALSE)))
}

#' Load all registry tables to list
#'
#' @return list
get_reg_ls <- function() {
  flatmultiread(c("E:/rawdata/709203/Eksterne data", "E:/rawdata/709203/Grunddata"))
}


#' Definition of relevant variables from DST tables
#'
#' @return
define_dst_vars <- function() {
  list(
    bef = c("PNR", "FAMILIE_ID"),
    faik = c("FAMILIE_ID", "FAMAEKVIVADISP_13", "FAMSOCIOGRUP_13"),
    ras = c("PNR", "SOC_STATUS_KODE"),
    uddf = c("PNR", "HFAUDD")
  )
}

#' Simple wrapper of dplyr::select
#'
#' @param data
#' @param vars
#'
#' @return
select_vars <- function(data, vars) {
  data |> dplyr::select({{ vars }})
}

#' Subset DST tables to only include relvant variables.
#'
#' @param ls List of all registry tables
#'
#' @return
get_dst_tables <- function(ls) {
  dst_vars <- define_dst_vars()
  dst_tbl <- toupper(names(dst_vars))
  ls.all <- purrr::map(seq_along(dst_vars), function(i) {
    ls.reg <- ls[grepl(paste0("^", dst_tbl[i]), names(ls))] |> purrr::map(select_vars, vars = dst_vars[[i]])
    names(ls.reg) <- paste0("y", stringr::str_extract(names(ls.reg), "[0-9]{4}"))
    ls.reg
  })
  names(ls.all) <- dst_tbl
  ls.all
}

#' Wrapper to generate string matching pattern for stringr::str_detect()
#'
#' @param data
#'
#' @return
match_str <- function(data) {
  paste0("[", paste0(data, collapse = ","), "]")
}

#' Wrapper to generate string matching pattern for grepl()
#'
#' @param data character vector
#'
#' @return
match_str_grepl <- function(data) {
  paste0("(", paste0(data, collapse = "|"), ")")
}

# get_dst_tables(ls)

#' Generate sequence of previous N length
#'
#' @param data numeric vector of length 1
#'
#' @return
#'
#' @examples
#' last5y(10)
#' last5y(c(10, 6, 3))
lastNy <- function(data, n = 5) {
  paste0("y", seq((data - n), data) - 1)
}

#' Filter PNR (cpr) across list elements
#'
#' @param data list of tibbles to pass through
#' @param index index number (PNR/CPR)
#'
#' @return tibble
filterCPRacross <- function(data, index) {
  data |>
    purrr::map(function(i) {
      i[i$PNR == index, ]
    }) |>
    purrr::list_rbind()
}

#' Summarise data from last 5 years prior to inclusion
#'
#' @param data list with
#' @param v.median variables to get median
#' @param v.latest variables to get latest
#' @param v.mean variables to get mean
#'
#' @return tibble
previousNyears <- function(data, data.clin, n.years = 5, v.median = NULL, v.latest = c("FAMSOCIOGRUP_13", "SOC_STATUS_KODE"), v.mean = c("FAMAEKVIVADISP_13")) {
  df.cpryear <- data.clin |> dplyr::transmute(
    CPR = PNR,
    year = as.numeric(format(as.Date(rdate), "%Y"))
  )

  seqs <- purrr::map(df.cpryear$year, lastNy, n = n.years)

  seq_along(seqs) |>
    purrr::map(function(i) {
      df <- data[c(seqs[[i]])] |>
        filterCPRacross(index = df.cpryear$CPR[i]) |>
        dplyr::group_by(PNR) |>
        dplyr::summarise(
          dplyr::across(tidyselect::any_of(v.latest), \(x) tail(x, n = 1), .names = "{.col}.latest"),
          dplyr::across(tidyselect::any_of(v.mean), \(x) mean(x, na.rm = TRUE), .names = "{.col}.{n.years}.mean"),
          dplyr::across(tidyselect::any_of(v.median), \(x) median(x, na.rm = TRUE), .names = "{.col}.{n.years}.median")
        )
    }) |>
    purrr::list_rbind()
}


#' Extract relevant and summarised data from BAF and FAIK
#'
#' @param data list of dst data tables
#'
#' @return tibble
#'
#' @examples
#' get_reg_ls() |>
#'   get_dst_tables() |>
#'   get_beffaikras()
get_beffaikras <- function(data, clin.data = get_clinical()) {
  data <- data[stringr::str_detect(match_str(c("BEF", "FAIK", "RAS")), names(data))] |> purrr::list_flatten()

  years <- stringr::str_extract(names(data), "y[0-9]{4}")

  years[duplicated(years)] |>
    purrr::map(grep, years) |>
    purrr::map(function(i) {
      data[c(i)] |> purrr::reduce(dplyr::full_join)
    }) |>
    purrr::set_names(years[duplicated(years)]) |>
    previousNyears(data.clin = clin.data)

  ## BEF
  ## # Befolkningsoversigt. Data skal bruges for at kunne udtrække husstandsindkomst.
  ## FAIK
  ## # Familieindkomst. Familie id skal flættes med ID fra BEF for hvert år for at tage hensyn til evt skifte i status.
  ## FAMAEKVIVADISP_13 er relevante variabel for ækvivaleret indkomst
  ## Der findes også familiesocioøkonomisk status. Gør som Sine. Be done with it!
}

#' Title
#'
#' @return
#' @export
#'
#' @examples
#' data <- read_edu_level()
#' data |> dplyr::count(ISCED)
read_edu_level <- function() {
  haven::read_dta("E:/Formater/SAS formater i Danmarks Statistik/STATA_datasaet/Disced/c_audd_level_l1l4_k.dta") |>
    dplyr::transmute(
      HFAUDD = start,
      ISCED = AUDD_LEVEL_L1L4_K
    )
}

haven::read_dta(
  "E:/Formater/SAS formater i Danmarks Statistik/STATA_datasaet/Disced/c_audd_level_l1l3_k.dta") |> 
  dplyr::count(AUDD_LEVEL_L1L3_K)

# ls <- get_reg_ls() |> get_dst_tables()

get_uddf <- function(ls) {
  ls |>
    purrr::pluck("UDDF") |>
    purrr::pluck(1) |>
    dplyr::mutate(HFAUDD = as.character(HFAUDD)) |>
    dplyr::left_join(read_edu_level()) |>
    dplyr::group_by(PNR) |>
    dplyr::summarise(ISCED = max(ISCED), .groups = "keep") |>
    dplyr::mutate(
      ISCED = as.numeric(ISCED),
      ISCED_lvl = dplyr::case_match(ISCED, 0:2 ~ "low",
        3:4 ~ "medium",
        5:9 ~ "high",
        .default = NA
      ),
      ISCED_bin = dplyr::case_match(ISCED, 0:3 ~ "low",
                                    4:9 ~ "high",
                                    .default = NA
      )
    )
}


#' Collects all relevant variables from DST tables
#'
#' @param ls list of all registry tables
#' @param df.clin clinical data set
#'
#' @return tibble
#' @examples
#' get_reg_ls() |> get_dst(df.clin = get_clinical())
get_dst <- function(ls, df.clin) {
  ls_dst <- get_dst_tables(ls)
  ls_dst |>
    ## BEFxFAIKxRAS
    get_beffaikras(clin.data = df.clin) |>
    dplyr::full_join(
      ## UDDF
      get_uddf(ls_dst)
    )
}

#' Eases pipe renaming of columns
#'
#' @param data tibble, list or other object, for which names() makes sense. see ?setNames
#' @param prefix prefix to add
#' @param exclude names not to modify
#' @param new.names character vector of all new names
#'
#' @return object of same class as data
#'
#' @examples
#' set_colnames(data = mtcars, prefix = "WOW", exclude = "mpg")
set_colnames <- function(data, new.names = NULL, prefix = NULL, exclude = c("PNR", "CPR"), prefix.sep = "_") {
  if (is.null(new.names)) {
    nms <- names(data)
  } else {
    nms <- new.names
  }

  if (is.null(prefix)) {
    nms.mod <- nms
  } else {
    nms.mod <- paste(prefix, nms, sep = prefix.sep)
  }

  setNames(
    object = data,
    nm = dplyr::if_else(stringr::str_detect(match_str(exclude), nms),
      nms,
      nms.mod
    )
  )
}


#' Store of variable names for data sub-setting
#'
#' @return list
#' @examples
#' define_variables()
define_variables <- function() {
  list(
    talos = c(
      "age",
      "reg_female",
      "nihss_0",
      "reg_trombolyse",
      "reg_trombektomi",
      # "rtreat",
      "pase_0",
      "reg_alone",
      "reg_bmi",
      # "reg_hojde",
      # "reg_vaegt_alt",
      "reg_smoker",
      "reg_more_alc",
      "reg_hyperten",
      "reg_diabetes",
      "reg_tidl_tci",
      "reg_atriefli",
      "reg_ami",
      "reg_perifer_arteriel"
    ),
    clin = c(
      "age",
      "reg_female",
      "nihss_0",
      "reg_trombolyse",
      "reg_trombektomi",
      # "rtreat",
      "rtreat_placebo"),
    lifestyle=c(
      "pase_0",
      "pase_4",
      "reg_alone",
      "reg_bmi",
      # "reg_hojde",
      # "reg_vaegt_alt",
      "reg_smoker",
      "reg_more_alc",
      "reg_hyperten",
      "reg_diabetes",
      "reg_tidl_tci",
      "reg_atriefli",
      "reg_ami",
      "reg_perifer_arteriel"),
    lifestyle.events=c(
      "pase_0",
      "pase_4",
      "reg_alone",
      # "reg_bmi",
      # "reg_hojde",
      # "reg_vaegt_alt",
      "reg_smoker",
      "reg_more_alc",
      "reg_hyperten",
      "reg_diabetes",
      # "reg_tidl_tci",
      "reg_atriefli",
      # "reg_perifer_arteriel",
      "reg_ami"),
    lifestyle.bmi=c(
      "reg_bmi"
    ),
    ses=c(
      # "soc_status",
      # "soc_status_work",
      "soc_status_nowork",
      # "fam_indk",
      "fam_indk_hl",
      # "fam_indk_high",
      # "fam_indk_low",
      # "edu_level",
      # "edu_high",
      # "edu_low",
      "edu_level_hl"
    ),
    assess.events = c(
      "who_4",
      "mdi_4",
      "mfi_gen_4",
      "mrs_4_above1",
      "time",
      "status",
      "event.include"
    ),
    assess.events.pre = c(
      "who_4",
      "mdi_4",
      "mfi_gen_4",
      "mrs_0_above0",
      "time",
      "status",
      "event.include"
    ),
    assess.pred = c(
      "who_0",
      "mrs_0_above0"
    ),
    extra = c(
      "soc_status",
      "pase_0",
      "pase_4" ,
      "event"
    ),
    cpr = "PNR",
    ssri=c("rtreat",
           "time.all",
           "status.all",
           "rdate",
           "enddate"),
    quartiles = c(
      "pase_0_quartile",
      "pase_4_quartile"
    )
  )
}

#' Get var names in vector from group names. Possibility to keep all vars for as log as possible. Can be supplied to `gtsummary` functions
#'
#' @param groups vector of group names. See names(define_variables()) for options
#'
#' @return
#' @export
#'
#' @examples
get_var_vec <- function(v.groups){
  define_variables()[{{ v.groups }}] |> purrr::list_c()
}

#' SUbsets dataset based on variable group names as defined
#'
#' @param vector character vector of category names
#'
#' @return character vector
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |> get_vars(c("universal", "events"))
get_vars <- function(data, vars.groups, vars.vec=NULL) {
  data |> dplyr::select(tidyselect::any_of(c("pase_0","pase_4",".imp",".id",get_var_vec(vars.groups),vars.vec)))
}



#' Collect all relevant data for the events analysis data set
#'
#' @param ls ls of tibbles
#'
#' @return tibble
#' 
#' @examples
#' ls <- targets::tar_read(list_filtered)
#' 
collectall <- function(ls) {
  purrr::pluck(ls, "clinical") |>
    dplyr::left_join(purrr::pluck(ls, "all_events") |> set_colnames(prefix = "event"), by = c("PNR" = "CPR")) |>
    dplyr::left_join(purrr::pluck(ls, "dst") |> set_colnames(prefix = "dst"), by = "PNR")
}


## Formatting for analysis

#' Function to cut and group PASE data
#'
#' @param data data set including pase_0 and _4
#'
#' @return tibble
pase_cutter <- function(data, pase.rev = TRUE, drop.pase = FALSE, drop.nas=FALSE) {
  data.classes <- class(data)
  if ("mids" %in% data.classes) {
    data <- data |> mice::complete(action = "long", include = TRUE)
  }

  data <- data |>
    dplyr::mutate(dplyr::across(.cols = c("pase_0", "pase_4"), \(i) {
      cut(x = i, breaks = quantile(pase_0, na.rm = TRUE), labels = 1:4, include.lowest = TRUE)
    }, .names = "{.col}_quartile")) |>
    dplyr::mutate(pase_change = factor(dplyr::case_when(
      pase_0_quartile == 1 & pase_4_quartile == 1 ~ "Persistently low",
      pase_0_quartile %in% 2:4 &
        pase_4_quartile %in% 2:4 ~ "Persistently high",
      pase_0_quartile == 1 &
        pase_4_quartile %in% 2:4 ~ "Increase",
      pase_0_quartile %in% 2:4 &
        pase_4_quartile == 1 ~ "Decrease"
    ), ordered = FALSE),
    pase_change=factor(pase_change,levels=c("Persistently low", "Increase", "Decrease", "Persistently high")))

  
  if (drop.pase) {
    data <- data |> dplyr::select(-tidyselect::all_of(c("pase_0_quartile", "pase_4_quartile", "pase_0", "pase_4")))
  }


  if (pase.rev) {
    data <- data |> dplyr::mutate(
      pase_change = factor(pase_change, levels = c("Persistently high", "Decrease", "Increase", "Persistently low"))
    )
  }

  if (drop.nas) {
    data <- data |>
      dplyr::filter(!is.na(pase_change))
  }
  
  
  if ("mids" %in% data.classes) {
    data |> mice::as.mids()
  } else {
    data
  }
}

# as.Date(data$event_date.event)
define_status_time <- function(data) {
  data |> dplyr::mutate(dplyr::across(c("rdate", "enddate", "event_date.event"), ~ as.Date(.x)),
    time = difftime(dplyr::if_else(is.na(event_date.event), date_cutter(), event_date.event), enddate) |> lubridate::time_length("years"),
    time.all = difftime(dplyr::if_else(is.na(event_date.event), date_cutter(), event_date.event), rdate) |> lubridate::time_length("years"),
    status = as.integer(!is.na(event_event.type)),
    time = dplyr::if_else(time > censor_cutter(), censor_cutter(), time),
    time.all = dplyr::if_else(time.all > censor_cutter(), censor_cutter(), time.all),
    status = dplyr::if_else(time > censor_cutter(), FALSE, status),
    status.all = dplyr::if_else(time.all > censor_cutter(), FALSE, status),
    # status= dplyr::if_else(status,1,0),
    event.include = time > 0
  )
}

# as.integer(c(TRUE,FALSE))

#' Grouping soc status
#'
#' @param data tibble
#'
#' @return tibble
group_soc_status <- function(data) {
  data |>
    dplyr::mutate(
      soc_status = factor(dplyr::case_when(
        soc_status < 200 ~ "work",
        soc_status == 200 ~ "off",
        # only ~4 in the data set off work
        soc_status >= 200 ~ "outside"
      )),
      soc_status_work = soc_status == "work",
      soc_status_nowork = !soc_status_work
    )
}

#' Correction of character NA
#'
#' @param data tibble
#' @param char.missing character vector of entries to consider as NA
#'
#' @return tibble
correct_na <- function(data, char.missing = "NA") {
  data |> dplyr::mutate(dplyr::across(dplyr::where(is.character), ~ dplyr::na_if(.x, char.missing)))
}

#' Formatting the complete data set
#'
#' @param data the merged raw data set
#'
#' @return tibble
#' @examples
#' ds <- targets::tar_read(df_all_data) |>
#'   data_formatting() |>
#'   subset_df("mdi") 
#'  ds |> skimr::skim()
#'  ds |> View()
data_formatting <- function(data) {
  to_logical <- grep(match_str_grepl(c("missings", "incompletes")), names(data))

  suppressWarnings(
    data |>
      correct_na() |>
      dplyr::mutate(dplyr::across(all_of(to_logical), ~ .x == "TRUE")) |>
      dplyr::mutate(
        # This uses the work-corrected score
        # pase_0 = dplyr::if_else(pase_score_missings_w_0 | is.na(talos_pase10_0), NA, pase_score_sum_w_0),
        # pase_4 = dplyr::if_else(pase_score_missings_w_4 | is.na(talos_pase10_4), NA, pase_score_sum_w_4),
        # Below is the plain PASE scor used according to the manual used with TALOS
        pase_0 = dplyr::if_else(pase_score_missings_0, NA,pase_score_sum_0),
        pase_4 = dplyr::if_else(pase_score_missings_4, NA,pase_score_sum_4),
        who_0 = as.numeric(talos_who07_0),
        who_4 = as.numeric(talos_who07_4),
        mrs_0 = factor(substr(talos_mrs01_0, 1, 1), ordered = FALSE),
        mrs_0_above0 = (as.numeric(mrs_0) - 1) > 0,
        mrs_4 = factor(substr(talos_mrs01_4, 1, 1), ordered = FALSE),
        mrs_4_above1 = (as.numeric(mrs_4) - 1) > 1,
        mdi_4 = as.numeric(talos_mdi12_4),
        mfi_gen_4 = as.numeric(talos_mfi_gen_4),
        nihss_0 = as.numeric(talos_nihss16_0),
        soc_status = dst_SOC_STATUS_KODE.latest,
        fam_indk = cut(dst_FAMAEKVIVADISP_13.5.mean,
          breaks = quantile(dst_FAMAEKVIVADISP_13.5.mean, probs = seq(0, 1, 1 / 3), na.rm = TRUE),
          ordered_results = FALSE,
          labels = c("low", "medium", "high"),
          include.lowest = TRUE
        ),
        fam_indk_bin = cut(dst_FAMAEKVIVADISP_13.5.mean,
                       breaks = quantile(dst_FAMAEKVIVADISP_13.5.mean, probs = seq(0, 1, 1 / 2), na.rm = TRUE),
                       ordered_results = FALSE,
                       labels = c("low", "high"),
                       include.lowest = TRUE
        ),
        fam_indk_hl=forcats::fct_rev(fam_indk),
        fam_indk_high = dplyr::if_else(fam_indk_bin=="high",TRUE,FALSE),
        fam_indk_low = dplyr::if_else(fam_indk_bin=="low",TRUE,FALSE),
        edu_level = factor(dst_ISCED_lvl, ordered = FALSE, levels = c("low", "medium", "high")),
        edu_high = dplyr::if_else(dst_ISCED_bin=="high",TRUE,FALSE),
        edu_low = dplyr::if_else(dst_ISCED_lvl=="low",TRUE,FALSE),
        edu_level_hl=forcats::fct_rev(edu_level),
        rtreat_placebo=rtreat=="Placebo",
        event = factor(event_event.type)
      ) |>
      group_soc_status() |>
      define_status_time() |> 
      pragmatic_imputation(vec = c("reg_trombolyse", "reg_trombektomi","reg_any_perf"))
  )
}


#' Title
#'
#' @param data
#'
#' @return
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |>
#'   events_ready() |>
#'   View()
events_ready <- function(data,v.groups=c("clin","lifestyle.events","ses", "assess.events"),vars.vec=NULL) {
  out <- data |>
    get_vars(vars.groups = v.groups,vars.vec=vars.vec) 
  
  
  if ("event.include" %in% names(out)){
    out <- out |>
    dplyr::filter(event.include) |>
    # dplyr::filter(!is.na(pase_0),!is.na(pase_4))|>
    dplyr::select(-tidyselect::all_of("event.include"))# |> 
    # labelling_data()
  }
  
  return(out)
}


#' Title
#'
#' @param data
#'
#' @return
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted)
#' data |> events_ready() |>
#'   View()
talos_ready <- function(data,v.groups=c("talos","ssri")) {
  data |>
    get_vars(vars.groups = v.groups) |>
    dplyr::mutate(inc_time=lubridate::time_length(difftime(enddate,rdate),"years"),
                  rtreat=factor(rtreat,labels=c("Active","Placebo"))) |> 
    dplyr::select(-rdate,-enddate) |> 
    dplyr::rename(status="status.all",
                  time="time.all") |> 
    labelling_data()
}


#' A good-enough approximation of the defined analysis-population in the TALOS protocol
#' 
#' In practice this includes patients in the study for more than ~35 days
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
talos_analysis_pop_filter <- function(data){
  data |>
    (\(.x){
      .x |> split(.x$rtreat)
    })() |> 
    purrr::map2(c(268,284),\(.x,.y){
      .x |> 
        dplyr::slice_max(inc_time,n = .y,with_ties = FALSE)
    }) |> 
    dplyr::bind_rows()
}


talos_analysis_pop_filter_imp <- function(data){
  data |>
    mice::complete(action = "long", include = TRUE) |> 
    (\(.x){
      split(.x,.x$.imp)
    })() |> 
    purrr::map(talos_analysis_pop_filter) |> 
    dplyr::bind_rows() |> 
    mice::as.mids()
}

events_table <- function(data,by){
      list(
        "Overall"=data,
        split(data,data[by])
        ) |> 
    purrr::list_flatten() |> 
    purrr::imap(\(.x,.i){
      .x |> 
        dplyr::summarise(
          group=.i,
          py = sum(time),
          events = sum(status),
          events_pr_100 = 100 * events / py
        )
    }) |> 
    dplyr::bind_rows() |> 
    setNames(c("group","Patient Years","Events","Events pr 100 patient years")) |> 
    tidyr::pivot_longer(-group) |> 
    tidyr::pivot_wider(names_from = group,values_from = value)|>
    gt::gt() |>
    gt::fmt_number(columns = -1, decimals = 1)
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
talos_imp <- function(data){
  fun_impute(data=data,outcome.vars = c("time","status"),
             ignore = c("inc_time","rtreat"),
             pragmatic.reg = FALSE,
             pase.mod = FALSE)
}

#' Title
#'
#' @param data 
#' @param v.groups 
#'
#' @return
#' @export
#'
#' @examples
events_ready_small <- function(data,v.groups=c("clin","lifestyle.events","ses", "assess.events")) {
  data |>
    get_vars(vars.groups = v.groups) |>
    dplyr::filter(event.include) |>
    dplyr::filter(!is.na(pase_0),!is.na(pase_4))|>
    dplyr::select(-tidyselect::all_of("event.include"))# |> 
  # labelling_data()
}

#' Title
#'
#' @param date
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted) 
#' 
#' data |>
#'   prediction_ready() |>
#'   View()
prediction_ready <- function(data,
                             var.grps=c("clin","lifestyle","ses", "assess.pred")) {
  data |>
    get_vars(var.grps)|>
    dplyr::filter(!is.na(pase_0),!is.na(pase_4))#|> 
    # labelling_data()
}

## Data inspection and exploration
##
##
#' Title
#'
#' @param data
#' @param subdf
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_all_data) |>
#'   subset_df() |>
#'   View()
subset_df <- function(data, subdf = "pase") {
  data[grepl(paste0("^(", paste("PNR", subdf, paste0("talos_", subdf), sep = "|"), ")"), names(data))]
}

#' Imputation as a function, includes "pragmatic imputation"
#'
#' @param data 
#' @param outcome.vars 
#' @param ignore 
#' @param pragmatic.reg 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted) |> events_ready() 
#' data |>  labelling_data() |> fun_impute()
fun_impute <- function(data, outcome.vars = c("status", "time"), ignore = NULL,pragmatic.reg=TRUE,pase.mod=FALSE) {
  # data |> mice::md.pattern()
  if (pragmatic.reg){
  data <- data |> pragmatic_imputation()
  }
  
  ## Excluding entries with missing outcome measures
  data <- data |> 
    dplyr::filter(!dplyr::if_any(tidyselect::all_of(c(outcome.vars)), ~ is.na(.x)))
  
  init <- data |>
    mice::mice(maxit = 0)

  meth <- init$method
  meth[ignore] <- ""

  pred <- init$predictorMatrix
  pred[, c(outcome.vars)] <- 0

  # data_out <- data |> mice::futuremice(
  #   pred = pred,
  #   method = meth,
  #   print = FALSE,
  #   parallelseed = 8123,
  #   use.logical = FALSE,
  #   maxit = 20,
  #   m = 10
  # )
  
  data_out <- data |> mice::mice(
    pred = pred,
    method = meth,
    print = FALSE,
    seed = 8123,
    maxit = 20,
    m = 10
  )
  
  if (pase.mod){
    data_out <- data_out |> pase_cutter_mids()
  }

  data_out
  
  # lattice::densityplot(imp_data)
  # Regarding EVENTS
  #
  # On inspection/eye-balling densityplots looks reasonable with the current settings
  #
}

#' Function to cut PASE in mids object
#'
#' @param data mids object
#'
#' @return mids object
#' @export
#'
pase_cutter_mids <- function(data){
data |>
  mice::complete(action = "long", include = TRUE) |> 
  pase_cutter(drop.pase = FALSE, drop.nas = TRUE)|> 
  mice::as.mids()
  }


#' Completes events data set, option to impute
#'
#' @param data
#' @param impute
#'
#' @return mids or tibble
#' @examples
#' targets::tar_read(df_all_data_formatted) |> events_dataset(impute=FALSE)
#' targets::tar_read(df_all_data_formatted) |> 
events_dataset <- function(data, impute = TRUE, uv=FALSE) {
  # data <- data |> events_ready()
  if (impute) {
    data |> events_ready()|>
      dplyr::select(-tidyselect::any_of("reg_bmi"))|>
      fun_impute(ignore = c("pase_0","pase_4"),pase.mod = FALSE) 
  } else if (uv){
    data |> events_ready()|>
      pase_cutter(drop.pase = TRUE,drop.nas = TRUE)
  } else if ("mids" %in% class(data)){
    data
  } else {
    data |> events_ready()|>
      dplyr::select(-tidyselect::any_of("reg_bmi")) |>
      pase_cutter(drop.pase = TRUE,drop.nas = TRUE)
  }
}

#' Title
#'
#' @param data 
#' @param all.vars 
#' @param outcome.var 
#' @param use.strata 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |> events_ready() |>  subset_df("pase")
#' data <- targets::tar_read(df_event_data) |> events_dataset(impute = FALSE) 

cox_regression <- function(data, all.vars = TRUE,outcome.var="pase_change",use.strata=TRUE, include_formula=FALSE) {
  if ("mids" %in% class(data)) {
    nms <- names(data$data)
  } else {
    nms <- names(data)
    data <- data |> 
      labelling_data()
    # BMI meassure is excluded from non-imputed dataset
    # data <- data |> dplyr::select(-tidyselect::all_of(c("reg_bmi")))
  }
  
  vars <- nms[!nms %in% c("time", "status", outcome.var)]

  form.prefix <- "survival::Surv(time, status) ~"
  
  if (use.strata) {
    reg.form <- glue::glue("{form.prefix} strata({outcome.var})")
  } else {
    reg.form <- glue::glue("{form.prefix} {outcome.var}")
    }

  if (all.vars) reg.form <- paste0(reg.form, " + ", paste(vars, collapse = " + "))
  
  require(survival)
  out <- with(data, survival::coxph(
    as.formula(reg.form)
  ))
  
  
  if (isTRUE(include_formula)){
    out$call$formula <- as.formula(reg.form)
  }
  
  out
}

#' Creates UV cox models for all variables in data set (but time and status)
#'
#' @param data 
#' @param include 
#'
#' @return
#' @export
#'
#' @examples
splitdf4uvcox <- function(data,include=c("pase_change","time","status")){
  names(data)[!names(data)%in%include] |> 
    purrr::map(\(.x){
      data |> dplyr::select(tidyselect::all_of(c(.x,include))) |> 
        cox_regression(outcome.var = .x,use.strata=FALSE,all.vars = FALSE)
               
      })
} 

#' Creates and prints UV Cox analyses
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_event_data) |> uv_cox_table()
uv_cox_table <- function(data){
  data |>
    events_dataset(impute = FALSE,uv = TRUE) |> 
    dplyr::select(pase_change,dplyr::everything()) |> 
  splitdf4uvcox(include=c("time","status")) |> 
  purrr::map(\(.x){
    .x |> tbl_regression_standard()
  }) |> 
    gtsummary::tbl_stack()
  }

tbl_merged_named <- function(data){
  data |> gtsummary::tbl_merge(tab_spanner = names(data))
}

tbl_stack_named <- function(data){
  data |> gtsummary::tbl_stack(group_header = names(data))
}

tbl_regression_standard <- function(data){
  # browser()
  data |> gtsummary::tbl_regression(exponentiate = TRUE, 
                                    add_estimate_to_reference_rows = TRUE,
                                    show_single_row = where(is.logical),
                                    statistics=list(gtsummary::all_continuous()~"[{conf.low};{conf.high}]",
                                                    gtsummary::all_categorical()~"[{conf.low}%;{conf.high}%]")
  ) |>
    fix_labels()
}

#' Wrapper to print summary table with extended info
#'
#' @param data formatted and subset data set
#' @param by.var stratify by
#'
#' @return
#' @examples
#' targets::tar_read(df_pred_data)|>print_table_summary(by="reg_female")
print_table_summary <- function(data, by.var = "pase_change") {
  data |> 
    labelling_data() |>
    # pase_cutter(drop.pase = TRUE) |>
    gtsummary::tbl_summary(
      missing = "ifany",
      by = tidyselect::all_of(by.var),
      value = list(where(is.logical) ~ TRUE)
    ) |>
    gtsummary::add_overall()
}

print_table_summary_explorer <- function(data, by.var = "pase_change") {
  data |> 
    labelling_data() |>
    # pase_cutter(drop.pase = TRUE) |>
    gtsummary::tbl_summary(
      missing = "ifany",
      by = tidyselect::all_of(by.var),
      value = list(where(is.logical) ~ TRUE),
      type = list(gtsummary::all_continuous() ~ "continuous2"),
      statistic = list(gtsummary::all_continuous() ~ c(
      #   # "{N_nonmiss} ({p_nonmiss}%)",
        "{median} ({p25}, {p75})",
      #   # "{min}, {max}",
        "{mean} ({sd})"#,
      #   # "{N_miss} ({p_miss}%)"
      )#,
      # gtsummary::all_categorical() ~ c(
      #   "{N_obs} ({p_nonmiss}%)"#,
      #   # "{N_miss} ({p_miss})"
      # )
      )
    ) |>
    gtsummary::add_overall() |>
    gtsummary::add_n() #|>
  # gtsummary::add_p()
}

#' Title
#'
#' @param data
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted)
#' targets::tar_read(df_all_data_formatted) |> events_tblone()
events_tblone <- function(data) {
  #browser()
  data |>
    # events_ready() |>
    pase_cutter(drop.pase = TRUE) |>
    dplyr::select(-tidyselect::all_of(c("status", "time"))) |>
    dplyr::filter(!is.na(pase_change)) |>
    print_table_summary()
}

#' Title
#'
#' @param data
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |> View()
#' targets::tar_read(df_pred_data)|>
#' dplyr::transmute(stRoke::quantile_cut(pase_0,4,group.names = 1:4),soc_status_work,fam_indk,edu_level) |> 
#' summary_tblone()
summary_tblone <- function(data,by=names(data)[1]) {
  # data <- targets::tar_read(df_pred_data)
  data |>
    labelling_data() |> 
    # prediction_ready() |>
    # dplyr::select(-reg_bmi) |>
    # pase_cutter(drop.pase = TRUE) |>
    # dplyr::filter(!is.na(pase_change)) |>
    # dplyr::mutate(pase_change=forcats::fct_rev(pase_change)) |> 
    print_table_summary(by.var = by)
}

#
#' Summaries of DST data for PASE quartiles at 0 and 4
#'
#' @param data 
#' @param vars 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_pred_data) |> sum_pase_tables()
sum_pase_tables <- function(data,vars=c("pase_0","pase_4")){
  vars |> lapply(function(.x){
    dplyr::tibble(stRoke::quantile_cut(data[[.x]],y=data[["pase_0"]],4,group.names = 1:4),
                  dplyr::select(data,soc_status_nowork,fam_indk_hl,edu_level_hl)) |> 
    summary_tblone()
    })
}


#' Creating a truthful stratified table for predictions
#'
#' @param data data frame
#'
#' @return list
#' @export
#'
#' @examples
#' targets::tar_read(df_pred_data) |> true_pred_sum_plot()
true_pred_sum_plot <- function(data){
  true_sum <-  data |>
    pase_cutter(drop.pase = TRUE) |>
    dplyr::mutate(pase_change=forcats::fct_rev(pase_change)) 
  
  list(#true_sum,
       true_sum |> (function(.x){
    split(.x,.x$pase_change %in% c("Persistently low","Increase"))
  })() |> 
    purrr::map(function(.y){
      .y |> 
        dplyr::mutate(pase_change=factor(pase_change))
      })) |> 
    purrr::list_flatten() |> 
    purrr::map(summary_tblone,by="pase_change") |>
    purrr::map(mask_micro_summary) |> 
    gtsummary::tbl_merge()
}

#' Get quick summary of missing vs non-missing for each given variable
#'
#' @param data data set
#' @param var variable to summarise over
#'
#' @return
#' @export
#'
#' @examples
#'  targets::tar_read(df_pred_data) |> dplyr::select(soc_status_work,fam_indk,edu_level) |> who_is_missing()
#'  targets::tar_read(df_pred_data) |> who_is_missing(var="reg_bmi")
who_is_missing <- function(data, var = "edu_level") {
  data |> 
    dplyr::mutate(log = factor(c("non-missing","missing")[is.na(data[[var]])+1])) |> 
    dplyr::select(log, tidyselect::everything(),-tidyselect::all_of(var)) |> 
    summary_tblone() #|> gtsummary::bold_p()
}

#' Title
#'
#' @param data
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted)
#' targets::tar_read(df_all_data_formatted) |> View()
#' targets::tar_read(df_pred_data) |> preds_tblone()
preds_tblone <- function(data) {
  data |>
    labelling_data() |> 
    # prediction_ready() |>
    # dplyr::select(-reg_bmi) |>
    pase_cutter(drop.pase = TRUE) |>
    dplyr::filter(!is.na(pase_change)) |>
    dplyr::mutate(pase_change=forcats::fct_rev(pase_change)) |> 
    print_table_summary()
}

#' Title
#'
#' @param data 
#' @param b.cols 
#'
#' @return
#' @export
#'
#' @examples
#' gt <- targets::tar_read(ls_pred_summary)[[1]]
#' gt |> add_var_groups_gt()
#' 
#' # For this to work, the function would need to handle labels and levels
#' gt <-  targets::tar_read(tbl_pred_summary)|> gtsummary::as_gt()
#' gt  |> add_var_groups_gt()
add_var_groups_gt <- function(gt){
  # gt <- sex_ls |> purrr::pluck(2) |> gtsummary::as_gt()
  # gt <- fix_labels(gt)
  cls <- class(gt)
  
  b.cols <- names(gt$`_data`)
  
  if (b.cols[[1]]!="variable"){
    # Flag to indicate if format is native gt or not. Simple assumption
    # class(gt) gt is not enough
    labels <- gt$`_data`[[1]]
    group.var <- names(gt$`_data`[[1]])
  } else {
    labels <- gt$`_data`[["label"]][gt$`_data`[["row_type"]]=="label"]
    group.var <- gt$`_data`[["variable"]]
  }
  
  
  
  groups <- matrix(ncol=length(labels)) |> 
    data.frame() |> 
    setNames(ifelse(labels=="","unknown_var",labels)) |> 
    tibble::as_tibble() |> groups_in_ds(labels = TRUE)
  
  group.labels <- names(groups) |> subset_named_labels(labels.raw = group_labels())

  labels.all  <-  group.labels |> purrr::imap(function(.x,.y){
    c(.x,groups[[.y]][["label"]])
  }) |> purrr::list_c()
  
  for (i in rev(seq_along(group.labels))){
    gt <- gt |> gt::tab_row_group(label=gt::md(glue::glue("*{group.labels[[i]]}*")),
                                  rows=which(group.var %in% groups[[names(group.labels)[[i]]]][["var"]]))
    
  }

  class(gt) <- cls
  gt
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(tbl_preds_lin_imp_reg)
#' data |> fix_labels()
fix_labels <- function(data){
  cls <- class(data)
  data[[1]][["variable"]][data[[1]][["row_type"]]=="label"] |> 
    subset_named_labels(var_labels()) |> 
    unname() -> data[[1]][["label"]][data[[1]][["row_type"]]=="label"]
  class(data) <- cls
  data
}


#' Title
#'
#' @param data 
#' @param b.cols 
#'
#' @return
#' @export
#'
#' @examples
#' tbl <- targets::tar_read(ls_pred_summary)[[1]]
add_var_groups_pre_calc <- function(data,b.cols){

  groups <- data |> groups_in_ds()
  
  group.labels <- names(groups) |> subset_named_labels(labels.raw = group_labels())
  
  t0 <- data.frame(matrix(ncol=length(b.cols))) |> 
    setNames(b.cols) |> 
    tibble::tibble()
  list(ext = group.labels |> purrr::imap(function(.x,.y){
    t0 |> dplyr::mutate(
      variable=.y,
      val_label=.x,
      row_type="group",
      label=.x
    )
  }) |> dplyr::bind_rows(),
  lvls = group.labels |> purrr::imap(function(.x,.y){
    c(.y,groups[[.y]][["var"]])
  }) |> purrr::list_c()
  )
}

#' Adds variable grouping and formatting to gtsummary tables
#'
#' @param tbl 
#'
#' @return
#' @export
#'
#' @examples
#' 
#' tbl <-  targets::tar_read(tbl_pred_summary)
#' targets::tar_read(tbl_pred_summary) |> add_var_groups()
#' 
add_var_groups <- function(tbl,
                           pre_ls=add_var_groups_pre_calc(tbl$inputs$data, 
                                                          names(tbl$table_body))){
  
  
  tbl |> gtsummary::modify_table_body(
    ~.x |> dplyr::bind_rows(pre_ls[["ext"]]) |> 
      dplyr::arrange(factor(variable,levels=pre_ls[["lvls"]]))
  ) |> 
    gtsummary::modify_table_styling(columns=label,
                                    rows= row_type%in%"level",text_format = "indent2") |>
    gtsummary::modify_table_styling(columns=label,rows= row_type%in%"label",text_format = "indent")|>
    gtsummary::modify_table_styling(columns=label,rows= row_type%in%"group",text_format = c("italic"))
}

#' Functionalised character vector of all labels
#'
#' @return
#' @export
#'
#' @examples
var_labels <- function(){
  c(
    age = "Age",
    reg_female = "Female sex",
    reg_bmi = "Body mass index",
    reg_smoker = "Current smoker",
    reg_alone = "Living alone",
    reg_more_alc = "High alcohol consumption",
    reg_hyperten = "Hypertension",
    reg_diabetes = "Diabetes",
    reg_atriefli = "Atrial fibrillation",
    reg_perifer_arteriel = "Peripheral arterial disease",
    reg_tidl_tci = "Previous TIA",
    reg_ami = "Previous MI",
    reg_trombolyse = "Treated with IVT",
    reg_trombektomi = "Treated with EVT",
    # reg_any_perf,
    rtreat = "Trial allocation",
    rtreat_placebo = "Placebo trial treatment",
    pase_0 = "Pre-stroke PASE score",
    pase_4 = "6 months post-stroke PASE score",
    pase_0_quartile = "Pre-stroke PASE score quartile",
    pase_4_quartile = "6 months post-stroke PASE score quartile",
    # pase_change,
    nihss_0 = "Admission NIHSS",
    # soc_status,
    soc_status_work = "Employed",
    soc_status_nowork = "Not employed",
    fam_indk = "Family income group",
    fam_indk_hl = "Lower family income",
    fam_indk_high = "Higher family income",
    fam_indk_low = "Lower family income",
    edu_level = "Educational level group",
    edu_level_hl = "Lower educational level",
    edu_high = "Higher educational level",
    edu_low = "Low educational level",
    who_4 = "WHO-5 score 6 months post-stroke",
    mdi_4 = "MDI score 6 months post-stroke",
    mrs_4_above1 = "mRS > 1 at 6 months post-stroke",
    mfi_gen_4 = "General fatigue (MFI domain) 6 months post-stroke",
    time = "Time",
    status = "Status",
    event.include = "Include event",
    who_0 = "Pre-stroke WHO-5 score",
    mrs_0_above0 = "Pre-stroke mRS > 0",
    pase_change = "PA change group"
  )
}


group_labels <- function(data){
  c("clin" = "Clinical data", 
    "lifestyle" = "Lifestyle and chronic diseases", 
    "ses" = "Socio-economic factors", 
    "assess.events" = "Assessments", 
    "assess.pred" = "Assessments", 
    "extra" = "extras")
}

rev_naming <- function(x){
  setNames(names(x),x)
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_pred_data)
groups_in_ds <- function(data, labels=FALSE){
  groups <- define_variables() |> purrr::imap(function(.x,.y){
    tibble::tibble(group=.y,var=.x)
  }) |> 
    dplyr::bind_rows()
  
  if (labels){
    matching <- subset_named_labels(names(data),
                        rev_naming(var_labels()))
  }else {
   matching <- names(data) 
  }
  groups[match(matching,groups[["var"]]),] |> 
    (\(.x){
      .x |> dplyr::mutate(group=factor(group,levels=unique(.x[["group"]])))
    })()  |> 
    cbind(
      tibble::tibble(
        label=labelling_data(data) |> labelled::var_label() |> purrr::list_c()
      )
    )|> 
    (\(.x){
    split(.x,.x[["group"]])
  })()
}


#' Subset labels
#'
#' @param data 
#' @param labels.raw 
#'
#' @return character vector
#' @export
#'
subset_named_labels <- function(data,labels.raw){
  labels.raw[match(data,names(labels.raw))]
}

#' Assign labels to data.frame or tibble
#'
#' @param data 
#' @param labels 
#'
#' @return
#' @export
#'
#' @examples
assign_labels <- function(data,labels){
  # data |> labelled::set_variable_labels(labels)
  
  labelled::var_label(data) <- labels

  data
}

#' Flexible labelling using labelled for nicer tables
#'
#' @param data data set
#'
#' @return
#' @export labelled data.frame/tibble
#'
#' @examples
#' data <- targets::tar_read(df_pred_data)
#' data <- data |> dplyr::mutate(test="test")
#' data |> labelling_data() |> labelled::var_label()
labelling_data <- function(data,label.list=var_labels()){

  labs <- subset_named_labels(names(data),label.list)
  labs[is.na(labs)] <- names(data)[is.na(labs)]
  
  data |> assign_labels(labels = labs)
}



#' Print regression table
#'
#' @param data cox regression ready data set
#'
#' @return gtsummary tbl_regression list object
#' @examples
#' data <- targets::tar_read(df_event_data)
#' targets::tar_read(df_events_mids) |> show_table_regression()
#' targets::tar_read(df_event_data) |> show_table_regression(use.mice=FALSE)
show_table_regression <- function(data, use.mice=FALSE, by.var="pase_change") {
  # browser()
  imp <- data |>
    events_dataset(impute = use.mice) 
  
  if ("mids" %in% class(imp)){
    imp <- imp |> pase_cutter_mids() |> 
      mice::complete(action = "long", include = TRUE) |> 
      dplyr::select(-tidyselect::any_of(c("pase_0","pase_4","pase_0_quartile","pase_4_quartile"))) |> 
      mice::as.mids()
  } else {
    imp <- imp  |> 
      dplyr::select(-tidyselect::any_of(c("pase_0","pase_4","pase_0_quartile","pase_4_quartile")))
  }
  imp |>  standard_multi_cox_table(by.var="pase_change")
}

standard_multi_cox_table <- function(data, by.var="pase_change",all.vars = TRUE){
  # browser()
  data |>
    cox_regression(all.vars = all.vars, use.strata = FALSE,outcome.var = by.var) |>
    tbl_regression_standard()
}


#' Splitting df to list by PA trajectory
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_pred_data) |> pred_ls_split()
pred_ls_split <- function(data, excluded.vars = "reg_bmi"){
  data |>
    pase_cutter(drop.pase = FALSE) |>
    dplyr::select(-tidyselect::all_of(c("pase_4","pase_0_quartile","pase_4_quartile"))) |> 
    dplyr::group_split(pase_split = pase_change %in% c("Increase", "Persistently low")) |>
    setNames(c("drop", "hop")) |>
    purrr::map2(.y = c("Decrease", "Increase"), .f = \(x, y){
      x |>
        dplyr::mutate(pase_bin = pase_change == y) |>
        dplyr::select(-tidyselect::all_of(c(excluded.vars, c("pase_change", "pase_split")))) |>
        na.omit()
    })
}

bin_original <- function(data,...){
  ## This will just follow the original pase_cutter binning
  data
}

#' Help developing new binning functions without using "browser()"
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted) |> bin_ready_boiler()
bin_ready_boiler <- function(data){
  data |>
    events_ready()|> 
    pase_cutter(drop.pase = FALSE)
}

#' Handles overall quantile change with dynamic definition of lowest and any change
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
bin_quantile <- function(data,quantiles=5,lowest=1,any=FALSE,...){
  
  highest_seq <- seq_len(quantiles)[-seq_len(lowest)]
  
  ## This will just follow the original pase_cutter binning
  df <- data|>
    dplyr::mutate(dplyr::across(.cols = c("pase_0", "pase_4"), \(i) {
      as.numeric(stRoke::quantile_cut(x=i,groups=quantiles,y=pase_0,na.rm = TRUE,group.names=seq_len(quantiles)))
    }, .names = "{.col}_quantile")) 
  
  if (any){
  out <- df |> 
    dplyr::mutate(pase_change = factor(dplyr::case_when(
      pase_0_quantile < pase_4_quantile ~ "up",
      pase_0_quantile > pase_4_quantile ~ "down",
      pase_0_quantile %in% seq_len(lowest) & pase_4_quantile %in% seq_len(lowest) ~ "low",
      pase_0_quantile %in% highest_seq &
        pase_4_quantile %in% highest_seq ~ "high"
    ), ordered = FALSE))
  } else {
    out <- df |> 
      dplyr::mutate(pase_change = factor(dplyr::case_when(
        pase_0_quantile %in% seq_len(lowest) &
          pase_4_quantile %in% highest_seq ~ "up",
        pase_0_quantile %in% highest_seq &
          pase_4_quantile %in% seq_len(lowest) ~ "down",
        pase_0_quantile %in% seq_len(lowest) & pase_4_quantile %in% seq_len(lowest) ~ "low",
        pase_0_quantile %in% highest_seq &
          pase_4_quantile %in% highest_seq ~ "high"
      ), ordered = FALSE))
  }
  
  
    out |> dplyr::mutate(pase_change=factor(pase_change,levels=c("high", "down", "up", "low")))
}

bin_percentage<-function(data,percentage,...){
  data |>
    dplyr::mutate(pase_0_cut = as.numeric(cut(pase_0,quantile(pase_0,probs = c(0,percentage/100,1),na.rm = TRUE),include.lowest = TRUE,labels = 1:2)),
                  pase_4_cut = as.numeric(cut(pase_4,quantile(pase_0,probs = c(0,percentage/100,1),na.rm = TRUE),include.lowest = TRUE,labels = 1:2)),
                  pase_change = dplyr::case_when(
                    pase_0_cut > pase_4_cut ~ "down",
                    pase_0_cut < pase_4_cut ~ "up",
                    pase_0_cut %in% 1 ~ "low",
                    pase_0_cut %in% 2 ~ "high"
                  ),
                  pase_change = factor(pase_change, levels= c("high", "down", "up", "low"))
    )
}

bin_anyupdown<-function(data,low.q=1,high.q=2:4,...){
data |>
    dplyr::mutate(pase_0_quartile=as.numeric(pase_0_quartile),
    pase_4_quartile=as.numeric(pase_4_quartile),
    pase_change = dplyr::case_when(
      pase_0_quartile > pase_4_quartile ~ "down",
      pase_0_quartile < pase_4_quartile ~ "up",
      pase_0_quartile %in% low.q ~ "low",
      pase_0_quartile %in% high.q ~ "high"
    ),
    pase_change = factor(pase_change, levels= c("high", "down", "up", "low"))
    )
}

bin_absupdown<-function(data,abs.bin,low.q=1,high.q=2:4,...){
  data |>
    dplyr::mutate(pase_dif=(pase_4-pase_0),
                  pase_change = dplyr::case_when(
                    pase_dif > (abs.bin) ~ "up",
                    pase_dif < -(abs.bin) | pase_0 == 0 ~ "down",
                    pase_0_quartile %in% low.q ~ "low",
                    pase_0_quartile %in% high.q ~ "high"
                  ),
                  pase_change = factor(pase_change, levels= c("high", "down", "up", "low")))
}

bin_relupdown<-function(data,rel.bin,low.q=1,high.q=2:4,...){
  data |>
    dplyr::mutate(pase_rel_dif=(pase_4-pase_0)/pase_0,
    pase_change = dplyr::case_when(
      pase_rel_dif > (rel.bin/100) ~ "up",
      pase_rel_dif < -(rel.bin/100) ~ "down",
      pase_0_quartile %in% low.q ~ "low",
      pase_0_quartile %in% high.q ~ "high"
    ),
    pase_change = factor(pase_change, levels= c("high", "down", "up", "low")))
}

#' Title
#'
#' @param data 
#' @param ... 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |> group_format(binning.fun=bin_clusterlcm)
bin_clusterlcm <- function(data,...){
  data|>
    #events_ready()|> 
    #pase_cutter(drop.pase = FALSE)|>
    #dplyr::filter(!is.na(pase_0),!is.na(pase_4)) |> 
    dplyr::select(-tidyselect::any_of(c("pase_0_quartile","pase_4_quartile","pase_split","pase_dif","pase_rel_dif")))|> 
    lcm_cluster(n.clusters = 4) |> 
    final_membership() |> 
    dplyr::mutate(pase_change=clust) |> 
    dplyr::select(-clust)
}

#' Multi grouping
#'
#' @param data 
#' @param binning.fun 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted)
#' targets::tar_read(df_all_data_formatted)  |> group_format(binning.fun = bin_relupdown,rel.bin=20)
#' targets::tar_read(df_all_data_formatted)  |> group_format(binning.fun = bin_original,var.excl=NULL)
#' targets::tar_read(df_all_data_formatted)  |> group_format(binning.fun = bin_clusterlcm,var.excl=c("pase_0_quartile","pase_4_quartile","pase_split","pase_dif","pase_rel_dif"))
group_format <- function(data,var.keep=NULL,binning.fun,var.excl=c("pase_0","pase_4","pase_0_quartile","pase_4_quartile","pase_split","pase_dif","pase_rel_dif"),...){
  out_all <- data |>
    events_ready()|> 
    pase_cutter(drop.pase = FALSE) |> 
    #bin_relupdown(rel.bin=20)
    binning.fun(...)|>
    dplyr::filter(!is.na(pase_0),!is.na(pase_4))
  
  exl.ndx <- !names(out_all) %in% var.excl
  if (all(exl.ndx)){
    out <- out_all
  } else {
    out <- out_all[exl.ndx]
  }
  
  return(out)
}

ls_pase_pred_clean<-function(data,excluded.vars){
data|>
    purrr::map(\(.x){
      .x |>
        dplyr::select(-tidyselect::any_of(c(excluded.vars, 
        c("pase_4","pase_0_quartile","pase_4_quartile","pase_change", "pase_split","pase_rel_dif")))) |>
        na.omit()
    })
}


pred_ls_split_alt <- function(data, excluded.vars = "reg_bmi",type="relupdown",rel.bin){
  #browser()
  df<-data |>
    pase_cutter(drop.pase = FALSE)

  if (type=="anyupdown"){
    ls_binned<-df |>
    bin_anyupdown()|>
    (\(.x){
      list(.x|>dplyr::mutate(pase_bin=pase_change=="down")|>dplyr::filter(!pase_0_quartile==1),
      .x|>dplyr::mutate(pase_bin=pase_change=="up")|>dplyr::filter(!pase_0_quartile==4))
    })()

    names<-c("Any quartile down","Any quartile up")
  }  else if (type=="relupdown"){
    ls_binned<-df |>
    bin_relupdown(rel.bin=rel.bin)|>
    (\(.x){
      list(.x|>dplyr::mutate(pase_bin=pase_change=="down"),
      .x|>dplyr::mutate(pase_bin=pase_change=="up"))
    })()

    names<-c(paste0("More than ",rel.bin,"% down"),paste0("More than ",rel.bin,"% up"))
  } 
  
    
  ls_binned  |>
    setNames(names)|>
    ls_pase_pred_clean(excluded.vars = excluded.vars)
}


#' Run regularisation steps for split data set
#'
#' @param data selected data set
#'
#' @return list
#'
#' @examples
#' data <- targets::tar_read(df_pred_data)
#' targets::tar_read(df_pred_data) |> pred_models(auto.l = TRUE,weighted = FALSE,rel.bin=20,split.type="relupdown")
pred_models <- function(data, split.type="pase_bin", excludes = "reg_bmi",rel.bin=50,...) {
  if (split.type=="pase_bin"){
  ls_split <- data |>
    pred_ls_split(excluded.vars = excludes) 
  } else {
    ls_split <- data |>
      pred_ls_split_alt(excluded.vars = excludes,type=split.type,rel.bin=rel.bin) 
  }
  
  ls <- ls_split|>
    purrr::map(\(.x) regularisation_steps(.x,...))

  class(ls) <- c("regular_list", class(ls))
  ls
}

cross_mean_median_exp_table <- function(data) {
  nms <- paste0("v", seq_len(ncol(data)))

  cross_calcs <- data |>
    as.data.frame() |>
    setNames(nms) |>
    dplyr::rowwise() |>
    dplyr::transmute(
      median = median(dplyr::c_across(tidyselect::all_of(nms))),
      medianOR = exp(median),
      mean = mean(dplyr::c_across(tidyselect::all_of(nms))),
      meanOR = exp(mean)
    )
  
  dplyr::tibble(names = rownames(data), cross_calcs) |> 
    dplyr::select(-tidyselect::all_of(c("mean","median")))
}


gather_coefs_step1 <- function(data) {
  data |>
    list3levelpluck(lvl1 = "model", lvl2 = "B") |>
    purrr::map(purrr::reduce, cbind)
}

gather_coefs <- function(data) {
  # imputed.list <- "mids_regular_list" %in% class(data)

  if ("mids_regular_list" %in% class(data)) {
    data_step1 <- data |>
      purrr::map(gather_coefs_step1) |>
      purrr::map(purrr::reduce, cbind)
  } else if ("regular_list" %in% class(data)) {
    data_step1 <- data |> gather_coefs_step1()
  } else {
    stop("The supplied list has to be class 'mids_regular_list' or 'regular_list'")
  }

  data_step1 |>
    purrr::map(cross_mean_median_exp_table) |>
    purrr::reduce(dplyr::full_join, by = "names", suffix = paste0("_", names(data)))
}


#' Merge and print model coefficients. Pools datafrom mids analyses.
#'
#' @param data list
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(ls_pred_models)
#' targets::tar_read(ls_pred_models) |> print_pred_coefs()
#' targets::tar_read(ls_pred_mids_reg) |> print_pred_coefs() |> add_var_groups_gt()
print_pred_coefs <- function(data) {
  # data <- targets::tar_read(ls_pred_mids_reg)
  # nms <- names(data)

  if ("mids_regular_list" %in% class(data)) {
    type.table <- "Pooled regularised models"
  } else if ("regular_list" %in% class(data)) {
    type.table <- "Single regularised model"
  } else {
    stop("The supplied list has to be class 'mids_regular_list' or 'regular_list'")
  }

  merged_tbl <- data |>
    gather_coefs() |> 
    ## Leaving out the intercept
    (function(.x) .x[-1,])()
  
  sel_mean_med <- colnames(merged_tbl)[!grepl(pattern = "OR",colnames(merged_tbl))][-1]
  sel_or <- colnames(merged_tbl)[grepl(pattern = "OR",colnames(merged_tbl))]
  
  news <- subset_named_labels(merged_tbl$names,var_labels())
  
  merged_tbl <- merged_tbl |> dplyr::mutate(names=dplyr::if_else(is.na(news),names,news))

  gt_merged_tbl <- merged_tbl|>
    gt::gt() |> 
    gt::fmt_number(decimals = 5)
  
  merged_tbl_log <- merged_tbl |> dplyr::mutate(dplyr::across(tidyselect::all_of(sel_or), ~.x!=1),
                              dplyr::across(tidyselect::all_of(sel_mean_med), ~.x!=0))
  
  for (j in colnames(merged_tbl)[-1]) {
    
    i <- merged_tbl_log[[j]]
    
    gt_merged_tbl <- gt_merged_tbl |> gt::tab_style(style = list(
      gt::cell_text(weight="bold")
    ),
    locations = gt::cells_body(
      columns=j,
      rows = i
    )
      )}

  for (i in names(data)) {
    gt_merged_tbl <- gt_merged_tbl |>
      gt::tab_spanner(label = i, columns = tidyselect::ends_with(i))
  }

  gt_merged_tbl |> gt::tab_spanner(
    label = type.table,
    columns = -1
  )
}

#' Calculates confusionMatrix from contingency tables. Pools if object class is .
#'
#' @param data
#'
#' @return list
#'
#' @examples
#' targets::tar_read(ls_pred_mids_reg) |> multi_table_cfm()
#' targets::tar_read(ls_pred_models) |> multi_table_cfm()
multi_table_cfm <- function(data) {
  # data <- targets::tar_read(ls_pred_mids_reg)
  if ("mids_regular_list" %in% class(data)) {
    data <- data |> purrr::map(\(x){
      x |>
        # Test tables are plucked
        # purrr::map(\(y) y |> purrr::pluck("model") |> purrr::pluck("cMatTest"))|>
        list3levelpluck(lvl1 = "model", lvl2 = "cMatTest") |>
        # All tables are add together
        purrr::reduce(\(i, j) i + j)
    })
  } else if ("regular_list" %in% class(data)) {
    data <- data |> list3levelpluck(lvl1 = "model", lvl2 = "cMatTest")
  } else {
    stop("The supplied list has to be class 'mids_regular_list' or 'regular_list'")
  }

  data |>
    purrr::map(caret::confusionMatrix)
}

#' Collect and summarise auc meassures. Pools if "mids_regular_list" object
#'
#' @param data
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(ls_pred_mids_reg) |> multi_auc_summary()
#' targets::tar_read(ls_pred_models) |> multi_auc_summary()
multi_auc_summary <- function(data) {
  if ("mids_regular_list" %in% class(data)) {
    data_step <- data |> purrr::map(\(x){
      x |>
        # Test tables are plucked
        list3levelpluck(lvl1 = "model", lvl2 = "auc_test") |>
        # All tables are add together
        purrr::reduce(c)
    })
  } else if ("regular_list" %in% class(data)) {
    data_step <- data |>
      list3levelpluck(lvl1 = "model", lvl2 = "auc_test") |>
      purrr::map(c)
  } else {
    stop("The supplied list has to be class 'mids_regular_list' or 'regular_list'")
  }

  data_step |>
    purrr::map(summary)
}

#' Map and 2 level recursive purrr::pluck to ease regular_list subsetting
#'
#' @param data
#' @param lvl1
#' @param lvl2
#'
#' @return
#' @export
#'
#' @examples
list3levelpluck <- function(data, lvl1 = "model", lvl2 = "cMatTest") {
  data |> purrr::map(\(y) y |>
    purrr::pluck(lvl1) |>
    purrr::pluck(lvl2))
}

#' Plot performance curve from glmnet regularisation
#'
#' @param data list of cvs.glmnet objects
#'
#' @return ggplot list object
#' @export
#'
#' @examples
plot_roc_curve <- function(data, title.text) {
  ggplot2::ggplot() +
    purrr::map(data, function(i) {
      ggplot2::geom_step(data = i, ggplot2::aes(x = FPR, y = TPR))
    }) +
    ggplot2::coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
    ggplot2::geom_abline() +
    ggplot2::theme_bw() +
    ggplot2::ggtitle(title.text)
}

roc_gather_step <- function(x) {
  with(x, glmnet::roc.glmnet(cvs[[1]]$fit.preval, newy = y1)[match(bestL, lambdas)])
}

#' Title
#'
#' @param data
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(ls_pred_mids_reg) |> multi_roc_plot()
#' targets::tar_read(ls_pred_models) |> multi_roc_plot()
multi_roc_plot <- function(data) {
  if ("mids_regular_list" %in% class(data)) {
    data_step1 <- data |>
      purrr::map(purrr::map, roc_gather_step) |>
      purrr::map(purrr::list_flatten)
  } else if ("regular_list" %in% class(data)) {
    data_step1 <- data |> purrr::map(roc_gather_step)
  } else {
    stop("The supplied list has to be class 'mids_regular_list' or 'regular_list'")
  }

  data_step1 |>
    purrr::map2(.y = names(data), plot_roc_curve) |>
    patchwork::wrap_plots()
}

#' Title
#'
#' @param tuning.param
#' @param data
#'
#' @return
#' @export
#'
#' @examples
multi_tuning_gather <- function(tuning.param = "bestA", data) {
  if ("mids_regular_list" %in% class(data)) {
    data_step1 <- data |>
      purrr::map(purrr::map, \(x) x |> purrr::pluck(tuning.param)) |>
      purrr::map(purrr::reduce, c)
  } else if ("regular_list" %in% class(data)) {
    data_step1 <- data |>
      purrr::map(purrr::pluck, tuning.param)
  } else {
    stop("The supplied list has to be class 'mids_regular_list' or 'regular_list'")
  }

  data_step1 |> purrr::map(summary)
}


#' Tidied tuning summary call
#'
#' @param data
#'
#' @return
#'
#' @examples
#' targets::tar_read(ls_pred_mids_reg) |> tuning_summary()
#' targets::tar_read(ls_pred_models) |> tuning_summary()
tuning_summary <- function(data) {
  c(ALPHA = "bestA", LAMBDA = "bestL") |> purrr::map(\(x) x |> multi_tuning_gather(data = data))
}

#' Apply regularisation steps to MIDS object, output arranged by grouping
#'
#' @param data mids object from mice package
#'
#' @return list
#'
#' @examples
#' targets::tar_read(df_pred_mids) |> mids_regularisation()
#' ls <- targets::tar_read(df_pred_mids) |> mids_regularisation(weighted = FALSE,rel.bin=20,split.type="relupdown")
mids_regularisation <- function(data,...) {
  ls <- data |>
    mice::complete(action = "long") |>
    dplyr::group_split(.imp) |>
    purrr::modify(\(x){
      x |> dplyr::select(-tidyselect::all_of(c(".imp", ".id")))
    }) |>
    purrr::map(\(.x)pred_models(.x,auto.l=TRUE,excludes = NULL,...))

  nms <- ls |>
    purrr::map(names) |>
    unique() |>
    purrr::reduce(c)

  # As a consequence of the above code each "set" of analyses are together.
  # Here the same group analyses are subset and grouped
  ls_n <- purrr::map(nms, function(i) {
    ls |> purrr::map(purrr::pluck, i)
  }) |>
    setNames(nms)

  # Special class is applied to ease future handling
  class(ls_n) <- c("mids_regular_list", class(ls_n))
  ls_n
}

#' A collection of all the summary functions to be applied to list of
#' pred_models() output
#'
#' @param data list of data
#'
#' @return list
#' @export
#'
multi_summary <- function(data){
  list( "coefTable" = print_pred_coefs(data) |> 
          gt::fmt_number(n_sigfig = 3) |> 
          fix_labels() #|> add_var_groups_gt()
        ,
        "confusionMatrices" = multi_table_cfm(data),
        "summaryAUC" = multi_auc_summary(data),
        "rocPlots" = multi_roc_plot(data),
        "tuningSummaries" = tuning_summary(data))
}

# funs <-list(
#     "coefTable" = print_pred_coefs,
#     "confusionMatrices" = multi_table_cfm,
#     "summaryAUC" = multi_auc_summary,
#     "rocPlots" = multi_roc_plot,
#     "tuningSummaries" = tuning_summary
# )

# multi_summary <- plyr::each(
#   "coefTable" = print_pred_coefs,
#   "confusionMatrices" = multi_table_cfm,
#   "summaryAUC" = multi_auc_summary,
#   "rocPlots" = multi_roc_plot,
#   "tuningSummaries" = tuning_summary
# )

#' Subset multiple elements from list
#'
#' @param data list
#' @param indices numeric or character vector
#'
#' @return list
#' @examples
#' targets::tar_read(ls_pred_summary)$confusionMatrices |> purrr::map(list_subset)
list_subset <- function(data,indices=c("overall","byClass")){
  data[indices]
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(ls_pred_summary) |> print_model_resutls()
print_model_resutls <- function(data){
  par(mfrow=c(1,2))
  
  list(data$coefTable,
       invisible(data$confusionMatrices |> purrr::map(\(x) x |> purrr::pluck("table"))),
       data$confusionMatrices |> purrr::map(list_subset),
       data$summaryAUC,
       data$tuningSummaries
       )
}


#' Classic logistic regression on prediction covariables
#'
#' @param data data frame
#'
#' @return
#' @export
#'
#' @examples gtsummary list elemnt
#' targets::tar_read(df_pred_data) |> pred_log_reg()
#' targets::tar_read(df_pred_data) |> pred_ls_split()
pred_log_reg <- function(data){
  data |> pred_ls_split() |> 
    purrr::map(\(x) {
      gtsummary::tbl_regression(glm(pase_bin~.,family = binomial,data = x),
                                exponentiate= TRUE)#|> 
        # gtsummary::bold_p()
    }
    ) |> (\(x){gtsummary::tbl_merge(tbls = x, 
                                    tab_spanner = names(x))})() }

#' Small wrapper to format CI with square brackets
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
square_ci <- function(data){
  gsub("(-?\\d*\\.?\\d*)(, )(-?\\d*\\.?\\d*)",
       "\\[\\1; \\3\\]",data)}

#' Apply CI formatting across gtsummary table including merged tables
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
fix_ci <- function(data){
  data |> gtsummary::modify_table_body(~ .x |> 
                                         dplyr::mutate(dplyr::across(dplyr::starts_with("ci"),
                                                       function(.y){square_ci(.y)})))
}

#' Classic linear regression on 6 months PASE score. Uni and multi.
#'
#' @param data data frame
#'
#' @return
#' @export
#'
#' @examples gtsummary list elemnt
#' data <- targets::tar_read(df_pred_data)
#' data <- targets::tar_read(df_pred_mids)
#' targets::tar_read(df_pred_data) |> pred_lin_reg()
#' targets::tar_read(df_pred_mids) |> pred_lin_reg()
pred_lin_reg <- function(data){
  
  # list("tbl_regression-str:ref_row_text"="Reference") |> 
  #   gtsummary::set_gtsummary_theme()
  
  if ("mids" %in% class(data)){
    cols <- names(data$data)
    
  } else {
    cols <- names(data)
    data <- data |> 
      labelling_data()
  }
  
  vars <- cols[cols!="pase_4"]
  
  formula_pase <- paste("pase_4",paste(vars,collapse = "+"),sep="~" )
  
  # multi <- with(data=data,lm(pase_4~.)) |>
  #   gtsummary::tbl_regression(add_estimate_to_reference_rows = TRUE)|> 
  #   gtsummary::bold_p() |> gtsummary::add_n()
  
  if (!"mids" %in% class(data)){
    ls <- list("Univariable"=data |> 
                 gtsummary::tbl_uvregression(method=lm, show_single_row = dplyr::where(is.logical),
                                             y=pase_4,
                                             add_estimate_to_reference_rows = TRUE,pvalue_fun = NULL)#|> gtsummary::bold_p()
               ,
               "Multivariable (no BMI)"=lm(pase_4~.,data=dplyr::select(data,-reg_bmi)) |>
                 gtsummary::tbl_regression(add_estimate_to_reference_rows = TRUE, show_single_row = dplyr::where(is.logical))|> 
                 # gtsummary::bold_p() |> 
                 gtsummary::add_n()
               ,
               "Multivariable (ALL)"= lm(pase_4~.,data=data) |>
                 gtsummary::tbl_regression(add_estimate_to_reference_rows = TRUE, show_single_row = dplyr::where(is.logical))|> 
                 # gtsummary::bold_p() |> 
                 gtsummary::add_n()
    )  
    
  } else {
    ls <- list("Multivariable (ALL)"= suppressWarnings(mice::lm.mids(pase_4~.,data=data) |>
                 gtsummary::tbl_regression(add_estimate_to_reference_rows = TRUE, show_single_row = dplyr::where(is.logical))|> 
                 # gtsummary::bold_p() |> 
                   gtsummary::add_n()))
  }
  
  ls |> (\(x){gtsummary::tbl_merge(tbls = x, 
                                   tab_spanner = names(x))})() |> 
    fix_ci()
}


#' Simple standard plot
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_event_data) |> events_dataset(impute = FALSE)|> cox_regression() |> plot_survival()
plot_survival <- function(data){
  data |>
    ggsurvfit::survfit2() |>
    ggsurvfit::ggsurvfit(linetype_aes = TRUE, size = 0.8) +
    ggsurvfit::add_confidence_interval() +
    ggsurvfit::add_risktable(
      risktable_stats = c("n.risk", "cum.event"),
      stats_label = list(cum.event = "Cumulative Observed Events",
                         n.risk = "Number at Risk"),
      theme =
        list(
          ggsurvfit::theme_risktable_default(axis.text.y.size = 11,
                                  plot.title.size = 11),
          ggplot2::theme(plot.title = ggplot2::element_text(face = "bold"))
        )
    ) +
    ggplot2::scale_y_continuous(
      limits = c(0, 1),
      labels = scales::percent,
      expand = c(0.01, 0)
    ) +
    ggplot2::scale_x_continuous(breaks = 0:9, expand = c(0.02, 0))
}


#' Smooth tidy survfit object
#'
#' @param data survfit object
#'
#' @return tibble
#' @export
#'
#' @examples
#' ls <- targets::tar_read(df_event_data) |> events_dataset(impute = FALSE)|> cox_regression()
#' data <- ls |> ggsurvfit::survfit2(robust=TRUE) |> 
#'   ggsurvfit::tidy_survfit(type="survival") |> 
#'   dplyr::group_split(strata) |> purrr::pluck(1)
#'   
#' ls |> ggsurvfit::survfit2(robust=TRUE) |> 
#'   ggsurvfit::tidy_survfit(type="survival") |> 
#'   dplyr::group_split(strata) |> 
#'   purrr::map(smooth_col)
smooth_col <- function(data, force_mono=TRUE){
  smoothed <- lapply(c("estimate","conf.high","conf.low"),function(i){
    # stats::predict(cobs::cobs(x = data$time,
    #                           y = data[i],
    #                           constraint = "decrease",
    #                           nknots=4,
    #                           pointwise = rbind(c(0,min(data$time),1)),
    #                           degree = 2,)) |> 
    #   tibble::as_tibble() |> dplyr::select(fit) |> 
    stats::predict(mgcv::gam(data=data,formula = as.formula(glue::glue("{i}~s(time,bs='cs')")))) |>
       tibble::as_tibble()|>
      setNames(glue::glue("{i}_smooth"))
    }) |> purrr::list_cbind()
 
  if (force_mono){
    ## Forcing starting point to be 1
    smoothed[1,1] <- 1
    
    smoothed[1] <- force_decrease(smoothed[1]) ## Only monotonize the esitimate
  }
  
  dplyr::tibble(data,
                smoothed)
   
}

#' Forces the direction
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
force_decrease <- function(data){
  data |> purrr::imap(function(.x,.n){
    s <- c()
    for (i in seq_along(.x)){
      if (i == 1) {
        s[1] <- .x[1]
      } else {
        if (.x[i]>s[i-1]){
          s[i] <- s[i-1]
        }  else {
          s[i] <- .x[i]
        }
      }
    }
    s
  }) |> 
    dplyr::bind_cols()
}


#' Prepare cox regression for smooth survival plot
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |> events_dataset(impute = FALSE)|> cox_regression() |> smooth_cox_data()
smooth_cox_data <- function(data){
  data |>
    ggsurvfit::survfit2(robust=TRUE) |> 
    ggsurvfit::tidy_survfit(type="survival") |> 
    dplyr::group_split(strata) |> 
    purrr::map(smooth_col) |> 
    purrr::list_rbind()
}

#' Plot smooth survival plot
#'
#' @param data df from cox regression
#'
#' @return ggplot list object
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_event_data) |> events_dataset(impute = FALSE) |> cox_regression(include_formula=TRUE)
#' data <- targets::tar_read(df_events_mids) |> cox_regression()
#' data |> plot_survival_smooth()
plot_survival_smooth <- function(data,line.w=1.2){
  if ("mira" %in% class(data)) stop("Only plots non-imputed survival data")
  
  n.level <- length(data$xlevels[[1]])
  
  # data |>  
  #   ggsurvfit::survfit2() |> 
  #   ggsurvfit::tidy_survfit()
  
  if (data |> ggsurvfit::survfit2() |> purrr::pluck("n") |> length() ==1 ){
    ds <- data |>  
      ggsurvfit::survfit2() |> 
      ggsurvfit::tidy_survfit() 
    p <- ds |> 
      ggplot2::ggplot(ggplot2::aes(x=time, y=estimate))+
      ggplot2::geom_smooth(se=TRUE, method="loess", formula = "y~x", linewidth=line.w, color="grey10")
    # Added auto max for y axis removed again to ensure same y axis
    # max_y <- max(ds$conf.high)
    
  } else {
    ds <- data |> 
      smooth_cox_data() 
    p <- ds |> 
      ggplot2::ggplot()+
      ggplot2::geom_line(ggplot2::aes(x=time, y=estimate_smooth, color=strata, linetype=strata), linewidth=line.w)+
      ggplot2::geom_ribbon(ggplot2::aes(x=time, ymin=conf.low_smooth,ymax=conf.high_smooth, fill=strata), alpha=.2)
    # Added auto max for y axis removed again to ensure same y axis
    # max_y <- max(ds$conf.high_smooth)
    
  }
  
  if (n.level==4){
    colors <- viridisLite::turbo(n=n.level,direction = 1)[c(1,3,2,4)]
  } else {
    colors <- viridisLite::turbo(n=n.level,direction = 1)
  }
  
  p+
    ggplot2::scale_y_continuous(limits = c(0,1.02),
      breaks = seq(0,1,.25),
      labels = scales::percent,
      expand = c(0.01, 0)
    ) +
    ggplot2::scale_x_continuous(breaks = 0:9, expand = c(0.02, 0))+
    ggplot2::scale_fill_manual(values=colors)+
    ggplot2::scale_color_manual(values=colors)+
    ggplot2::theme_minimal()+
    ggplot2::theme(axis.title.x = ggplot2::element_blank(),
                   axis.title.y = ggplot2::element_blank(),
                   # axis.text = ggplot2::element_blank(),
                   # legend.position = "none",
                   panel.grid.minor.y = ggplot2::element_blank(),
                   panel.grid.major.y = ggplot2::element_line(color="grey45",linewidth = line.w/2))
  
}

# viridisLite::turbo(n=4,direction = 1)[c(1,3,2,4)]

cluster_rank <- function(data){
  data |> cox_regression(outcome.var = "clust",use.strata = TRUE) |> ggsurvfit::survfit2(robust=TRUE) |> 
    ggsurvfit::tidy_survfit(type="survival") |> 
    dplyr::group_split(strata) |> 
    purrr::map(\(x){
      min(x[["estimate"]])
    }) |> purrr::list_c() |> rank() |> rev()
}

cox_relevel <- function(data){
  data |> dplyr::mutate(clust=factor(clust,levels=cluster_rank(data)))
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_all_data_formatted)
complete_preds_data <- function(data){
  data |> 
    # events_ready() |>
    fun_impute(ignore = c("pase_0","pase_4"),pase.mod = FALSE) |> 
    mice::complete() |> 
    dplyr::filter((!is.na(pase_0)&!is.na(pase_4)))
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |> add_kamila_cluster()
#' data_kam <- targets::tar_read(df_all_data_formatted) |> add_kamila_cluster() 
#' data_kam |>print_table_summary(by.var = "kam_grp")
#' data_kam  |> cox_regression(outcome.var="kam_grp")|> plot_survival_smooth()
#' data_kam  |> cox_regression(outcome.var="kam_grp",use.strata = FALSE)|> gtsummary::tbl_regression(exponentiate = TRUE, add_estimate_to_reference_rows = TRUE) |> gtsummary::bold_p()
#' targets::tar_read(df_events_complete) |> kamila_cluster(n.clusters=3)
kamila_cluster <- function(data, n.clusters=3,include.out=FALSE){
  # An index number could be added to later join pack. Of input a complete data set from imputation and pooling??
  data_orig <- data
  
  
  if (!include.out){
    data <- data |> 
      dplyr::select(-tidyselect::one_of(c("time","status")))
  }
  
  
  catInd <- data |> lapply(\(x) is.character(x)|is.logical(x)) |> purrr::list_c()
  conInd <- data |> lapply(\(x) is.numeric(x)|is.integer(x)) |> purrr::list_c()
  
  catVars <- data[,catInd] 
  catVars <- catVars |> lapply(factor) |> dplyr::bind_cols() |> as.data.frame()
  conVars <- data[,conInd] |> scale()|> as.data.frame()
  
  if (is.null(n.clusters)){
    out <- kamila::kamila(conVar = conVars, catFactor = catVars, numClust = 2:7, numInit = 10,
                             calcNumClust = "ps"
    )
  }else {
  out <- kamila::kamila(conVar = conVars, catFactor = catVars, numClust = n.clusters, numInit = 10)
  }
  
  ls <- list("out"=out,"data_orig"=data_orig)
  
  class(ls) <- c("kamila_cluster",class(ls))
  
  ls
  
}


#' VarSelLCM wrapper
#'
#' @param data complete dataset with no missings
#' @param n.clusters number of clusters (if length 1, n is fixed, in n>1 given clusters are tested)
#' @param include.out flag to include outcome variables or not
#' @param memb.out output data frame with final membership or not (then outputs standard model output) 
#'
#' @return list with VarSelLCM output and original dataset with cluster appended
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_events_complete)
#' data |> lcm_cluster()
lcm_cluster <- function(data, n.clusters=3, include.out=FALSE, var.sel=FALSE,out.vars=c("time", "status")){
  data_orig <- data
  
  if (!include.out){
    data <- data |> 
      dplyr::select(-tidyselect::any_of(out.vars))
  }
  
 
  set.seed(5432)
  
  out <- data |> 
    dplyr::mutate(dplyr::across(where(is.logical)|where(is.character),~factor(.x))) |> 
    as.data.frame() |> 
    VarSelLCM::VarSelCluster(
    gvals=n.clusters,
    crit.varsel="BIC",
    vbleSelec = var.sel,
    nbcores = round(parallel::detectCores()*.8)
    )
  
  ls <- list("out"=out,"data_orig"=data_orig)
  
  class(ls) <- c("lcm_cluster",class(ls))
  
  ls
  
}



#' Kmeans clustering
#'
#' @param data 
#' @param n.clusters 
#' @param include.out 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_events_complete)
#' data |> kmeans_cluster()
#' data |> kmeans_cluster(n.clusters=3)
kmeans_cluster <- function(data, n.clusters=3, include.out=FALSE, memb.out=TRUE){
  data_orig <- data
  
  if (!include.out){
    data <- data |> 
      dplyr::select(!tidyselect::one_of(c("time", "status")))
  }
  
  out <- data |>
    dplyr::mutate(dplyr::across(where(is.double),~scale(.x)),
      dplyr::across(where(is.logical)|where(is.character),~factor(.x)),
      dplyr::across(where(is.factor),~as.numeric(.x))) |> 
    stats::kmeans(
      centers=n.clusters
    )
  
  ls <- list("out"=out,"data_orig"=data_orig)
  
  class(ls) <- c("kmeans_cluster",class(ls))
  
  ls
  
}

#' dbscan clustering
#'
#' @param data 
#' @param n.clusters 
#' @param include.out 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_events_complete)
#' data |> dbscan_cluster()
#' data |> dbscan_cluster(n.clusters=3)
dbscan_cluster <- function(data, n.clusters=3, include.out=FALSE, memb.out=TRUE){
  data_orig <- data
  
  if (!include.out){
    data <- data |> 
      dplyr::select(!tidyselect::one_of(c("time", "status")))
  }
  
  data <- data |> na.omit() |> dplyr::mutate(rtreat=rtreat!="Placebo",
                                             dplyr::across(dplyr::everything(), as.numeric)) 
  
  
  ## This plot indicates that eps should be set around 60, but at this value everything is one cluster.
  dbscan::kNNdistplot(data,k = 5)
  
  ## Performing hierachical clustering, it is clear, that the algorithm is not able to seperate clusters.
  hds <- dbscan::hdbscan(data,minPts = 5)
  
  plot(hds,show_flat = TRUE)
  
  ## Clustering with set eps value and minPts
  ds <- dbscan::dbscan(data,eps = 25,minPts = 2)
  
  ds[["cluster"]]
  
  ## dbscan is not an interesting approach, apparently
  
  # 
  # 
  # 
  # 
  # out <- data |>
  #   dplyr::mutate(dplyr::across(where(is.double),~scale(.x)),
  #                 dplyr::across(where(is.logical)|where(is.character),~factor(.x)),
  #                 dplyr::across(where(is.factor),~as.numeric(.x))) |> 
  #   stats::kmeans(
  #     centers=n.clusters
  #   )
  # 
  # ls <- list("out"=out,"data_orig"=data_orig)
  # 
  # class(ls) <- c("kmeans_cluster",class(ls))
  # 
  # ls
  
}

#' Title
#'
#' @param ls 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_events_complete)
#' ls <- data |> lcm_cluster()
#' ls |> final_membership()
final_membership <- function(ls){
  cls <- class(ls)
  if ("kamila_cluster" %in% cls) {
    
      tibble::tibble(clust=factor(ls$out$finalMemb),
                     ls$data_orig)
    
  } else if ("lcm_cluster" %in% cls) {
    
    tibble::tibble(clust=factor(ls$out@partitions@zMAP),
                   ls$data_orig)
    
  } else if ("kmeans_cluster" %in% cls) {
    
    tibble::tibble(clust=factor(ls$out$cluster),
                   ls$data_orig)
    
  } else stop("Class not recognised")
  
}


#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_events_complete)
#' data |> get_clusters()
#' data |> get_clusters(n.cl=2:7)
get_clusters <- function(data,n.cl=4,rm.out=TRUE){
  set.seed(1123)
  
  if (length(n.cl)>1){
    list(
      # "kmeans"=data |> kmeans_cluster(n.clusters = n.cl,include.out = !rm.out),
      "lcm"=   data |> lcm_cluster(n.clusters = n.cl,include.out = !rm.out),
      "kamila"=data |> kamila_cluster(n.clusters = n.cl,include.out = !rm.out)
    )
  } else {
    list(
      "kmeans"=data |> kmeans_cluster(n.clusters = n.cl,include.out = !rm.out),
      "lcm"=   data |> lcm_cluster(n.clusters = n.cl,include.out = !rm.out),
      "kamila"=data |> kamila_cluster(n.clusters = n.cl,include.out = !rm.out)
    )
  }
  
}

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(list_pred_clusters)
#' data |> final_clusters()
final_clusters <- function(data,new.levels=NULL){
  out <- data |> lapply(final_membership)  |> lapply(labelling_data)
  
  if (is.null(new.levels)){
    out
  } else {
    out  |> 
      purrr::map2(relevels,function(x,y){
        # x$clust <- factor(factor(x$clust,levels=y),labels=1:4)
        x$clust <- factor(x$clust,levels=y)
        x #|> 
        # dplyr::filter(clust %in% range(as.numeric(clust))) |> 
        # dplyr::mutate(clust=factor(clust))
      })
  }
  
  }



#' Title
#'
#' @param data 
#' @param by 
#'
#' @return
#' @export
#'
#' @examples
merged_summary_tbl <- function(data,by="clust"){
  data |>
    purrr::map(function(x){
      x |> 
        # print_table_summary(by.var = by)
        gtsummary::tbl_summary(by=by) #|>
        # gtsummary::add_p() |> gtsummary::bold_p()
    }
    ) |> (\(x){
      x |> gtsummary::tbl_merge(tab_spanner = names(x))
    })()
}

#' Easy cox regression tbl for uniform results
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
cox2tbl <- function(data,by="clust",all.vars=FALSE){
data|> cox_regression(outcome.var = by,use.strata = FALSE,all.vars = all.vars) |>
  gtsummary::tbl_regression(exponentiate =TRUE,
                            statistics=list(gtsummary::all_continuous()~"[{conf.low};{conf.high}]",
                                            gtsummary::all_categorical()~"[{conf.low}%;{conf.high}%]")) |> 
    gtsummary::bold_p()
}

#' Title
#'
#' @param data 
#' @param by 
#'
#' @return
#' @export
#'
#' @examples
merged_cox_reg_tbl <- function(data,by="clust"){
  data |>
    purrr::map(function(x){
      x |> cox2tbl(by=by)
    }
    ) |> (\(x){
      x |> gtsummary::tbl_merge(tab_spanner = names(x))
    })()
}

#' Title
#'
#' @param data 
#' @param by 
#'
#' @return
#' @export
#'
#' @examples
wrapped_surv_plot <- function(data,by="clust"){
data |>
  purrr::map(function(x){
    x |> cox_regression(outcome.var = by,use.strata = TRUE,all.vars = FALSE) |> 
      plot_survival_smooth()+ggplot2::labs(color="Cluster",fill="Cluster",linetype="Cluster")
  }
  ) |> (\(x){
    x |> patchwork::wrap_plots(ncol=1) + patchwork::plot_annotation(tag_levels = list(names(x)))
  })()
}


# Ranking by most events
relevel_by_rank <- function(data){
  ## Assigning clusters to each dataset
data <- targets::tar_read(list_pred_clusters) |> 
  final_clusters(new.levels = NULL)

## Calculating cox regressions and ranking by the final point on the survival plot
relevels <- data |>
  purrr::map(function(x){
    x |> cox_regression(outcome.var = "clust",use.strata = TRUE,all.vars = FALSE) |>  
      ggsurvfit::survfit2() |> 
      ggsurvfit::tidy_survfit() |> 
      (\(x){
        split(x,x[["strata"]]) |> 
          purrr::map(function(.y){
            .y[["estimate"]][nrow(.y)]
          })
      })() |> purrr::reduce(c) |> rank()
  }
  )

#3 Assigning the new, ranked levels
targets::tar_read(list_pred_clusters) |> 
  final_clusters(new.levels = relevels)
}

## TODO
## Verify definitions
## Do remaining documentation of functions
##
##
## How does elastic net work with imputed dataset?
## Functionalise to allow for imputed and non-imputed (both analyses) - in both cases with and without BMI - include department of inclusion to investigate reason of missing BMI data
##
## tidymodels does not allow pmm in mice. Thy're out!

## 

missing_stats <- function(index,data,glue.mask= "{N_miss} ({round(p_miss*100,dec)}%)",dec){
  index.var <- unique(data$table_body$variable)[index]
  
  table_body <- data$table_body |> 
    dplyr::filter(variable == index.var)
  
  if ("missing" %in% table_body$row_type){
    f1 <- grep("stat_\\d+",names(table_body))
    
    ## This approach 
    masks <- data$meta_data$df_stats |> 
      purrr::pluck(index)  |> 
      dplyr::filter(!duplicated(col_name)) |> 
      dplyr::arrange(col_name) |> 
      dplyr::mutate(mask=glue::glue(glue.mask))|> 
      dplyr::pull(mask)
    
    table_body[table_body$row_type=="missing",f1] <- masks |>
      as.matrix() |>
      t() |>
      tibble::as_tibble(.name_repair = "unique_quiet")
    
  } 
  table_body
}

missing_stats_steps <- function(body,ls,glue.mask,dec){
  seq_along(unique(body$variable)) |> 
    purrr::map(function(.x) {
      missing_stats(index=.x,data=ls,glue.mask = glue.mask,dec=dec)
    }) |> 
    dplyr::bind_rows()
}

add_missing_stats <- function(data,glue.mask= "{N_miss} ({round(p_miss*100,dec)}%)",dec=1){
  data |> 
    gtsummary::modify_table_body(
      ~ .x |> missing_stats_steps(ls=data,glue.mask = glue.mask,dec=dec)
    )
}

variable_masks <- function(index, data, cut.off,glue.mask= "<{n} (<{p}%)",dec=dec) {
  index.var <- unique(data$table_body$variable)[index]
  
  table_body <- data$table_body |> 
    dplyr::filter(variable == index.var)
  
  
  ## Filtering bu two different approaches
  if (table_body$var_type[1] %in% c("dichotomous", "categorical")) {
    if (any(grepl("^stat_[1-9]|[1-9]\\d",names(table_body)))){
    masked <- micro_n_masks(
      ## Handling nominal/binary
      
      body = table_body |> 
        dplyr::filter(row_type!="missing"),
      n.all = data$meta_data$df_stats |> 
        purrr::pluck(index) |> 
        dplyr::select(n),
      N.all=data$meta_data$df_stats |> 
        purrr::pluck(index) |> 
        dplyr::select(N),
      cut.off=cut.off,
      glue.mask = glue.mask
    )
  } else {
    masked <- table_body |> dplyr::filter(row_type!="missing")
  }
    
    
    if ("stat_0" %in% names(masked)){
      body <-  masked
      tb  <-  body |> dplyr::filter(row_type!="missing",!is.na(stat_0))
      f1  <-  grep("^stat_0", names(tb))
      
      meta.index <- data$meta_data$df_stats |> 
        purrr::pluck(index)
      
      if ("by" %in% names(meta.index)) {
        meta.index <- meta.index |> dplyr::filter(is.na(by))
        }
      
      masked <- masking(body=body,
                   tb=tb,
                   f1=f1,
            ns =  meta.index  |> 
              dplyr::select(n)|>
              dplyr::slice(seq_len(length(f1) * nrow(tb))) |>
              unlist(use.names = FALSE) |>
              matrix(ncol = length(f1), byrow = TRUE) |>  
              tibble::as_tibble(.name_repair = "unique_quiet"),
            Ns= meta.index |> 
              dplyr::select(N_obs) |> 
              dplyr::slice(1)|>
              unlist(use.names = FALSE),
            cut.off=cut.off,
            glue.mask=glue.mask,
            dec=dec)
    }
    
    out <- rbind(
      masked,
      table_body |> dplyr::filter(row_type=="missing"))
    
  } else {
    out <- table_body
  }
  
  if ("missing" %in% out$row_type) {
    ## Handling missings n is N_miss, and N is N_obs
    
    missings <- micro_n_masks(
      body = out |> 
        dplyr::filter(row_type == "missing"),
      n.all = data$meta_data$df_stats |> 
        purrr::pluck(index) |> 
        dplyr::select(N_miss),
      N.all=data$meta_data$df_stats |> 
        purrr::pluck(index) |> 
        dplyr::select(N_obs),
      cut.off=cut.off,
      glue.mask=glue.mask 
    )
    
    if ("stat_0" %in% names(missings)) {
      ## Handling overall column
      body  <-  missings
      tb <-  body |> dplyr::filter(row_type=="missing")
      f1 <-  grep("^stat_0", names(tb))
      
      masking(body=body,
              tb=tb,
              f1=f1,
              ns = data$meta_data$df_stats |> 
                purrr::pluck(index) |> 
                dplyr::filter(is.na(by)) |> 
                dplyr::select(N_miss) |> 
                dplyr::slice(1) |>
                unlist(use.names = FALSE) |>
                matrix(ncol = length(f1), byrow = TRUE) |>
                tibble::as_tibble(.name_repair = "unique_quiet"),
              Ns=data$meta_data$df_stats |> 
                purrr::pluck(index) |> 
                dplyr::filter(is.na(by)) |> 
                dplyr::select(N_obs) |> 
                dplyr::slice(1) |>
                unlist(use.names = FALSE),
              cut.off=cut.off*2,
              glue.mask=glue.mask,
              dec=dec)
      
    }
    
    out <- rbind(
      out |> 
        dplyr::filter(row_type != "missing"),
      missings)
    
    
    
  
    
  }
  
  out
}

micro_n_masks <- function(body, n.all, N.all, cut.off, glue.mask,dec=1) {
  if (nrow(body) > 1) {
    # First row removed in case of categorical
    # Possibly change to include in filter, to have whole df, or just rbind in the end
    tb <- body |> dplyr::filter(row_type!="label")
  } else { # last option is "dichotomous"
    tb <- body
  }

  ## Supports any number of stat columns (overkill!)
  f1 <-  grep("^stat_[1-9]|[1-9]\\d", names(tb))
  
  masking(body=body,
          tb=tb,
          f1 = f1,
          ns=n.all |>
            dplyr::slice(seq_len(length(f1) * nrow(tb))) |>
            unlist(use.names = FALSE) |>
            matrix(ncol = length(f1), byrow = TRUE) |> 
            tibble::as_tibble(.name_repair = "unique_quiet"),
          Ns=N.all |>
            dplyr::slice(seq_len(length(f1))) |>
            unlist(use.names = FALSE),
          cut.off=cut.off,
          glue.mask=glue.mask,
          dec=dec
          )
  
}


#' Title
#'
#' @param body full table body
#' @param tb filtered table body
#' @param f1 Subsets indexes of relevant columns
#' @param ns Relevant ns arranged in matrix
#' @param Ns All Ns 
#' @param cut.off 
#' @param glue.mask 
#' @param dec 
#'
#' @return
#' @export
masking <- function(body,
                    tb=NULL,
                    f1,
                    ns,
                    Ns,
                    cut.off=cut.off,
                    glue.mask=glue.mask,
                    dec=dec
                    ){
  
  if (is.null(tb)) tb <- body
  
  ## Logical matrix of relevant ns to mask
  f2 <- ns |>
    purrr::map(\(.x) .x %in% 1:(cut.off - 1)) |>
    dplyr::bind_cols() |>
    as.matrix()
  
  ## Number of cells with zero observations in each matrix row
  n0 <- apply(ns == 0, 1, sum)
  
  ## Number of cells with ns to mask in each matrix row
  ns.low <- apply(f2, 1, sum)
  
  
  if (all(ns.low==0)){
    out <- body
  } else {
    
    ls.rows <- seq_len(nrow(f2)) |>
      purrr::map(\(.i){
        
        ds <- tb[.i, ]
        # Creating masked matrix to record maskings
        masked <- matrix(FALSE, ncol = ncol(f2))
        
        # Only modify in case of low, handle one col matrices
        if (apply(f2, 1, any)[.i]) {
        n <- cut.off
        N <- Ns[f2[.i,]]
        p <- round(100 * n / N, dec)
        
        ds[f1[f2[.i,]]] <- glue::glue(glue.mask) |>
          as.matrix() |>
          t() |>
          tibble::as_tibble(.name_repair = "unique_quiet")
        masked[f2[.i,]] <- TRUE
        
        ## loop to add masks until satisfied
        while (sum(masked)==1 & # The case of only one low
               nrow(masked)>1 | # But ignored in case of ncol==1
               sum(ns[.i,][masked]) <= cut.off & 
               (sum(masked) + n0[.i]) < ncol(masked)) {
          # in the case that sum of smalls is below cutoff, another field is added.
          
          ranked <- apply(ns[.i,], 1, rank, ties.method = "first") |> t()
          
          ranked.i <- ranked == sum(masked) + n0[.i] + 1
          n <- ns[.i,][ranked.i] |> plyr::round_any(accuracy = cut.off, f = ceiling)
          N <- Ns[ranked.i]
          p <- round(100 * n / N, dec)
          
          ds[f1[ranked.i]] <- glue::glue(glue.mask)
          masked[ranked.i] <- TRUE
        }}
        
        list(
          ds = ds,
          masked = masked
        )
      })
    
    # The case for dichotomous
    if (nrow(tb) == 1) {
      out <- ls.rows |>
        purrr::map(purrr::pluck, "ds") |>
        dplyr::bind_rows()
      
    # Handling categorical data
    } else if (nrow(tb) > 1) {
      masks <- ls.rows |>
        purrr::map(purrr::pluck, "masked") |>
        purrr::reduce(rbind)
      
      out <- ls.rows |>
        purrr::map(purrr::pluck, "ds") |>
        dplyr::bind_rows()
      
      ## Indices by row
      col.i <- seq_len(nrow(masks)) |> purrr::map(\(.j){
        which(masks[.j,])
      })
      
      col.i.vec <- purrr::list_c(col.i) |> unique()
      
      if (purrr::compact(col.i) |> length() == 1){
        # As this is only the case with overall column
        # This should be reworked
        
        # This was the approach, to just select the first
        # which(purrr::map_lgl(col.i,is_empty))[1]
        
        # This will select the cell with the second lowest number
        col.i[[which(rank(ns,ties.method = "first")==2)]] <- c("")
      }
      
      out <- col.i |>
        purrr::map(\(.y){
          # length(.y)
          if (length(.y) > 0) {
            cols <- col.i.vec[!col.i.vec %in% .y]
            if (length(cols)==0){
              cols <- "" 
            } else {
              cols
            }
          } else {
            .y
          }
        }) |>
        purrr::imap(\(.y, .i){
          ds <- out[.i, ]
          if (length(.y) > 0 & all(.y!="")) {
            n <- ns[.i, .y] |> 
              purrr::map_dfr(plyr::round_any,accuracy = cut.off, f = ceiling)
            N <- Ns[.y]
            p <- round(100 * n / N, dec)
            if (n==0) glue.mask <- "{n}"
            
            ds[f1[.y]] <- glue::glue(glue.mask)|>
              as.matrix() |>
              t() |>
              tibble::as_tibble(.name_repair = "unique_quiet")
            ds
          } else {
            ds
          }
        }) |>
        dplyr::bind_rows()
      
      out <- rbind(
        body |>
          dplyr::filter(row_type=="label"),
        out
      )
    }
  }
  out
  
}


summary_masks <- function(body,ls,cut.off=5,dec=dec){
  seq_along(unique(body$variable)) |> 
    purrr::map(function(.x) {
      variable_masks(index=.x,data=ls,cut.off=cut.off,dec=dec)
    }) |> 
    dplyr::bind_rows()
}

mask_micro_summary <- function(data,micro.n=5){
  data |> 
    gtsummary::modify_table_body(
      ~ .x |> summary_masks(ls=data,cut.off = micro.n,dec=1)
    )
}

#' Title
#'
#' @param mask 
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' minimal_mask(mask=tibble::as_tibble(matrix(c(FALSE,FALSE,TRUE,TRUE),nrow=1),.name_repair="unique_quiet"),
#' data=tibble::as_tibble(matrix(c(20,12,8,4),nrow=1),.name_repair="unique_quiet"),micro.n=5)
minimal_mask <- function(mask,data,micro.n){
  ## Function assumes monotonous data
  ## Use with column_diffs() for risk tables
  
  ## FUnction works in practice, but example doesn't??
  
  out <- mask
  
  for (i in seq_len(nrow(mask))){
    if (i>1){
      for (j in seq_len(ncol(mask))){
        if(dplyr::pull(mask[i,j])){
          n <- max(which(which(!out[i,])<j))
          out[i,j] <- abs(data[i,n]-data[i,j]) %in% 1:(micro.n-1)
        }
      }
      }
    
  }
  
  out
}


#' Title
#'
#' @param data 
#' @param micro.n 
#' @param down 
#' @param col.sel 
#'
#' @return
#' @export
#'
#' @examples
#' data <- targets::tar_read(df_event_data) |>
#'   events_dataset(impute = FALSE) |>
#'   cox_regression(all.vars = FALSE) |>
#'   ggsurvfit::survfit2() |>
#'   ggsurvfit::tidy_survfit(times = c(0, 2, 4, 6, 8.5))|>
#'   tidyr::pivot_wider(id_cols = strata, names_from = time, values_from = cum.event)
#' data |> mask_micro_table(col.sel=-strata)
#' data.sel <- data |> dplyr::select(-strata)
mask_micro_table <- function(data, micro.n = 5, down = TRUE, col.sel) {
  data.sel <- data |>
    dplyr::select({{ col.sel }})
  
  ## List of the two selection matrices
  masked <- list(
    data.sel, # The actual data
    data.sel |>
      column_diffs(include.first = TRUE) # Row differences (incl first row)
  ) |>
    purrr::map(\(.y){ # For each element in the list, do colwise test
      .y |>
        purrr::map_dfr(\(.x) .x %in% 1:(micro.n-1))
    }) |>
    purrr::imap(\(.y,.i){ # apply minimal masking to last list object
      if (.i==2){
        .y |> minimal_mask(data = data.sel,micro.n=micro.n)
      }else {
        .y
      }
    }) |>
    purrr::reduce(`|`) |> # Combine matrices
    tibble::as_tibble() |> # To tibble
    purrr::map2(data.sel, \(.x, .y){ # Apply masking based on combined selection
      ifelse(.x, rounded_interval(.y, round = micro.n, down = down), .y)
    }) |>
    dplyr::bind_cols()
  
  ## Bind masked data to original columns
  dplyr::bind_cols(
    data |>
      dplyr::select(-{{ col.sel }}),
    masked|>
      # Converts new to character for uniform data
      dplyr::mutate(dplyr::across(dplyr::everything(), ~ as.character(.x)))
  ) |>
    dplyr::select(colnames(data)) # Order columns as original input data
}

rounded_interval <- function(data, round = 5, down = TRUE) {
  # Handle "ties"
  sub <- ifelse(down, -1, 1)
  data <- ifelse(data %% round == 0 & data != 0, data + 1, data)
  
  c(floor, ceiling) |>
    purrr::map(\(.x) {
      plyr::round_any(x = data, accuracy = round, f = .x)
    }) |>
    dplyr::bind_cols(.name_repair = "unique_quiet") |>
    setNames(c("l", "h")) |>
    dplyr::transmute(mask = glue::glue("{l}-{h}")) |>
    dplyr::pull(mask)
}

column_diffs <- function(data, prefix.pattern = NULL, include.first = TRUE, suffix.out = "_diff") {
  if (!is.null(prefix.pattern)) {
    data <- data |>
      dplyr::select(tidyselect::starts_with(prefix.pattern))
  }
  
  index <- seq_along(data)[-1]
  
  diff <- index |>
    purrr::map(\(.y){
      abs(data[.y - 1] - data[.y])
    }) |>
    dplyr::bind_cols() |>
    (\(.x) setNames(.x, paste0(names(.x), suffix.out)))()
  
  if (include.first) {
    out <- dplyr::bind_cols(data[1], diff)
  } else {
    out <- diff
  }
  out
}


collect_calibration<-function(data){
  data|>
  purrr::map(\(.y) {
    .y |>
      purrr::pluck("model") |>
      purrr::pluck("TestProb") |>
      dplyr::bind_rows()
      })
}

collect_calibration_mids<-function(data){
  data|>
  purrr::map(\(.z) {
    .z |>
      purrr::map(\(.y) {
        .y |>
          purrr::pluck("model") |>
          purrr::pluck("TestProb") |>
          dplyr::bind_rows()
      }) |> 
        dplyr::bind_rows()
        })
}

plot_calibration<-function(data,name) {
          predtools::calibration_plot(
            data = as.data.frame(
              dplyr::select(data, y, pred) |>
                dplyr::mutate(y = as.numeric(y) -
                                1)
            ),
            obs = "y",
            pred = "pred"#,
            # x_lim = c(0, 1),
            # y_lim = c(-.1, 1.1)
          ) |>
        purrr::pluck("calibration_plot") +
        ggplot2::labs(title = name)+
        ggplot2::scale_x_continuous(breaks=seq(0,1,.25),limits=c(0, 1))+
        ggplot2::scale_y_continuous(breaks=seq(0,1,.25),limits = c(-.1, 1.1))
}

print_calibration<-function(data,file){
  ggplot2::ggsave(filename =file,plot = data,device = "png",dpi = 600,width = 84,height = 150,units = "mm")
}

map_summary_results <- function(data){
  data |> lapply(\(.x){
    .x |>
      multi_summary() |>
      print_model_resutls()
  })
}


#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' list("Raw dec/inc quartiles" = targets::tar_read(ls_pred_models)) |> map_summary_results() |> minimal_performance_table()
minimal_performance_table <- function(data){
  data |>
    purrr::imap(\(.x, .i){
      df_perf <- lapply(.x[[3]], \(.y){
        Reduce(c, .y)
      }) |> dplyr::bind_rows()
      
      df_ns <- .x[[2]] |>
        purrr::map(\(.y) {
          tibble::tibble(N=sum(.y),
                         n=sum(.y[2,]))
        }) |>
        dplyr::bind_rows()
      
      
      df_auc <- .x[[4]] |>
        purrr::map(\(.y) .y[["Median"]]) |>
        purrr::reduce(c)
      
      ## Here should follow a block to add median alpha and lambda to the table overview
      
      df_tune <- purrr::map(.x[[5]], \(.y){
        .y |>
          purrr::map(\(.z){
            .z[["Median"]]
          }) |>
          purrr::reduce(c)
      }) |> dplyr::bind_cols()
      
      tibble::tibble(
        group = .i,
        model = names(.x[[3]]),
        df_ns,
        df_perf,
        "Median AUC" = df_auc,
        setNames(df_tune, paste("Median", names(df_tune)))
      )
    }) |>
    dplyr::bind_rows() |>
    dplyr::select(group, model, N,n,Sensitivity, Specificity, `Pos Pred Value`, `Neg Pred Value`, tidyselect::starts_with("Median"))
}

data_prep_mmrm <- function(data){
  data |>
    dplyr::mutate(id = factor(dplyr::row_number())) |>
    dplyr::select(id, dplyr::everything(), -reg_bmi) |>
    tidyr::pivot_longer(
      cols = dplyr::starts_with("pase_"),
      values_to = "pase", names_to = "time"
    ) |>
    dplyr::mutate(time = factor(time), ) |>
    dplyr::mutate(dplyr::across(dplyr::where(is.logical), \(.x) as.numeric(.x)))
}

split_sex <- function(data,var="reg_female"){
  data |> 
    (\(.x){
      split(dplyr::select(.x,-tidyselect::all_of(var)),.x[[var]]) |> 
        setNames(c("male","female"))
    })()
}

simple_multi_mmrm <- function(data,vars.out=c("pase", "time", "id")){
  vars <- names(data)[!names(data) %in% c("pase", "time", "id")]
  
  ## mmrm doesn't work too well with gtsummary as variable sorting is lost
  mmrm::mmrm(as.formula(paste0("pase~", paste(vars, collapse = "+"), "+us(time|id)")), data = data)
}

simple_uni_mmrm <- function(data,vars.out=c("pase", "time", "id")){
  vars <- names(data)[!names(data) %in% c("pase", "time", "id")]
  
  ## mmrm doesn't work too well with gtsummary as variable sorting is lost
  vars |> purrr::map(\(.x){
    mmrm::mmrm(as.formula(paste0("pase~", .x, "+us(time|id)")), data = data)
  })
  
}

mmrm_summary <- function(data){
  gtsummary::tbl_regression(data,
                            add_header_row = TRUE,
                            show_single_row = tidyselect::where(is.logical),
                            tidy_fun = broom.helpers::tidy_parameters
  ) |>
    fix_labels() |> 
    gtsummary::modify_table_styling(column = p.value,
                                    hide=TRUE)
}


df_mega_list <- function(
    strategies=c(
  "bin_original", 
  "bin_anyupdown", 
  "bin_relupdown", 
  "bin_absupdown",
  "bin_quantile", 
  "bin_clusterlcm",
  "bin_percentage")){
    purrr::map(strategies, \(.z){
      if (.z == "bin_original") {
        f <- eval(str2expression(.z))
        list(c(list(binning.fun = f), list(var.excl = c("pase_4", "pase_0_quartile", "pase_4_quartile", "pase_split", "pase_dif", "pase_rel_dif")), list(name = "original_change")))
      } else if (.z == "bin_anyupdown") {
        f <- eval(str2expression(.z))
        purrr::map(c("_q1", "_half"), \(.y){
          if (.y == "_q1") {
            out <- list(binning.fun = f, low.q = 1, high.q = 2:4)
          } else if (.y == "_half") {
            out <- list(binning.fun = f, low.q = 1:2, high.q = 3:4)
          }
          c(out, list(var.excl = c("pase_4", "pase_0_quartile", "pase_4_quartile", "pase_split", "pase_dif", "pase_rel_dif")), list(name = glue::glue("any_change{.y}")))
        }) # |> purrr::list_flatten()
      } else if (.z == "bin_relupdown") {
        f <- eval(str2expression(.z))
        purrr::map(c(20, 50, 75), \(.x){
          out <- purrr::map(c("_q1", "_half"), \(.y){
            if (.y == "_q1") {
              out <- list(binning.fun = f, rel.bin = .x, low.q = 1, high.q = 2:4)
            } else if (.y == "_half") {
              out <- list(binning.fun = f, rel.bin = .x, low.q = 1:2, high.q = 3:4)
            }
            c(out, list(var.excl = c("pase_4", "pase_0_quartile", "pase_4_quartile", "pase_split", "pase_dif", "pase_rel_dif")), list(name = glue::glue("rel_change_{.x}{.y}")))
          })
        }) |> purrr::list_flatten()
      } else if (.z == "bin_absupdown") {
        f <- eval(str2expression(.z))
        purrr::map(seq(40, 100, 20), \(.x){
          out <- purrr::map(c("_q1", "_half"), \(.y){
            if (.y == "_q1") {
              out <- list(binning.fun = f, abs.bin = .x, low.q = 1, high.q = 2:4)
            } else if (.y == "_half") {
              out <- list(binning.fun = f, abs.bin = .x, low.q = 1:2, high.q = 3:4)
            }
            c(out, list(var.excl = c("pase_4", "pase_0_quartile", "pase_4_quartile", "pase_split", "pase_dif", "pase_rel_dif")), list(name = glue::glue("abs_change_{.x}{.y}")))
          })
        }) |> purrr::list_flatten()
      } else if (.z == "bin_percentage") {
        f <- eval(str2expression(.z))
        purrr::map(c(5, 15), \(.x){
           out <- list(binning.fun = f, percentage = .x)
           c(out, list(var.excl = c("pase_4", "pase_0_quartile", "pase_4_quartile", "pase_split", "pase_dif", "pase_rel_dif","pase_0_cut","pase_4_cut")), list(name = glue::glue("lowest_percentage_{.x}")))
          })
      } else if (.z == "bin_quantile") {
        f <- eval(str2expression(.z))
        purrr::map(seq(6, 10, 2), \(.x){
          purrr::map(seq_len(floor(.x/2)), \(.y){
            purrr::map(c(FALSE,TRUE),\(.any){
              out <- list(binning.fun = f, quantiles = .x, lowest = .y, any = .any)
              
              c(out, 
                list(var.excl = c("pase_4", "pase_0_quartile", "pase_4_quartile", "pase_split", "pase_dif", "pase_rel_dif" 
                                  ,
                                  "pase_0_quantile", "pase_4_quantile"
                                  )), 
                list(name = glue::glue("quantile_change_{.x}_{.y}{ifelse(.any,'_ANY','')}"))
                )
              
            })
            
          })|> purrr::list_flatten()
        }) |> purrr::list_flatten()
      } else if (.z == "bin_clusterlcm") {
        f <- eval(str2expression(.z))
        list(c(
          list(binning.fun = f),
          list(var.excl=c("pase_0_quartile","pase_4_quartile","pase_split","pase_dif","pase_rel_dif")), #Returns all false, which crashes the function
          list(name = "cluster_lcm")
        ))
      }
    }) |> purrr::list_flatten()
}


multi_grouping_wrapper <- function(name, ...) {
  .f <- function(...) {
    data |>
      group_format(...)
  }
  
  list(.f(...)) |> setNames(name)
}

# ls <- df_mega_list()

# ls <- c(list(binning.fun = bin_quantile, quantiles = 5, lowest = 2, any=FALSE), list(var.excl = c("pase_4", "pase_0_quartile", "pase_4_quartile", "pase_split", "pase_dif", "pase_rel_dif")), list(name = glue::glue("quantile_change_5_2")))
# # , "pase_0_quantile", "pase_4_quantile"
# multi_grouping_df_list(data=targets::tar_read(df_all_data_formatted),list(ls)) |> purrr::pluck(1) |> View()

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(df_all_data_formatted) |> multi_grouping_df_list(args.list=df_mega_list())
multi_grouping_df_list <- function(data,args.list,remove=NULL){
  purrr::map(args.list, \(.x){
    if ("mids" %in% class(data)){
      out <- data |>
        mice::complete(action = "long", include = TRUE) |> 
        (\(.y){
          do.call(multi_grouping_wrapper, c(list(data = .y), .x))
        })() |> 
        purrr::map(\(.z){
          .z |>
            dplyr::select(-tidyselect::any_of(remove)) |> 
            mice::as.mids()
        }) 
    } else {
      out <- do.call(multi_grouping_wrapper, c(list(data = data), .x))|> 
        purrr::map(\(.z){
          .z |>
            dplyr::select(-tidyselect::any_of(remove))
        }) 
    }
  out
  }) |> purrr::list_flatten()
  }

#' Title
#'
#' @param data 
#'
#' @return
#' @export
#'
#' @examples
#' targets::tar_read(list_df_multi_grouping)[24] |> multi_results_list()
multi_results_list <- function(data){
  data |> 
    purrr::map(\(.x){
      list(summary=.x |> gtsummary::tbl_summary(by = pase_change),
           cox = .x |> 
           # dplyr::select(-pase_0,-pase_4,-pase_rel_dif) |>
           cox_regression(all.vars = TRUE, use.strata = FALSE, outcome.var = "pase_change") |>
             tbl_regression_standard(),
           plot=.x |>
             # dplyr::select(-pase_0,-pase_4,-pase_rel_dif) |>
             cox_regression() |>
             ggsurvfit::survfit2() |>
             ggsurvfit::ggsurvfit() + ggsurvfit::add_confidence_interval() +
             ggsurvfit::scale_ggsurvfit() +
             ggsurvfit::add_risktable()
      )
    })}

multi_cox_performance_test <- function(data,...){
  data |>
    purrr::imap(\(.x,.i){
      out <- .x |>
        cox_regression(...) |>
        performance::model_performance()
      tibble::tibble(Name=.i,out)
    })|> 
    dplyr::bind_rows() |> 
    dplyr::mutate(rank=rank(AIC,ties.method = "min"))
    # performance::compare_performance(rank = TRUE)
}


mids_model_aic <- function(data){
  sapply(data[["analyses"]],AIC) |> 
    median()
} 

pick_non_duplicated <- function(data,name,aic,i){
  data[[name]][!duplicated(data[[aic]])][i]
}