4131 lines
117 KiB
R
Executable file
4131 lines
117 KiB
R
Executable file
# 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]
|
|
}
|