PAaSO/defence.R

source(here::here("1 PA Decline/dst import.R"))

source(here::here("fun/functions.R"))

################################################################################
############ Study 1
################################################################################

# pak::pak("agdamsbo/project.aid")

############ MMRM

ls <- project.aid::docx2list(path = here::here("1 PA Decline/Fra DDV/241002/mmrm.docx"))

df_mmrm_coef <- setNames(ls[[1]][c(1, 4, 5)], c("var", "beta", "ci"))

df_mmrm_coef <- df_mmrm_coef[-nrow(df_mmrm_coef), ]

df_mmrm_coef_clean <- apply(df_mmrm_coef, 2, \(.x) {
  .x[.x == "—"] <- "0, 0"
  .x
}) |>
  as.data.frame() |>
  tidyr::separate(
    col = "ci", into = c("lo", "hi"), sep = ", ", convert = TRUE
  ) |>
  apply(2, \(.x) {
    .x[.x == "0, 0"] <- "0"
    .x[.x == ""] <- NA
    .x
  }) |>
  as.data.frame() |>
  tibble::as_tibble() |>
  dplyr::mutate(dplyr::across(tidyselect::all_of(c("beta", "lo", "hi")), as.numeric))

## Manual cleaning
df_mmrm_coef_clean$var <- c(
  "Age", "Female sex", "Admission NIHSS", "Treated with IVT",
  "Treated with EVT", "Placebo trial treatment", "Living alone",
  "Current smoker", "High alcohol consumption", "Hypertension",
  "Diabetes", "Previous TIA", "Atrial fibrillation", "Previous MI",
  "Peripheral arterial disease", "Not employed", "Lower family income",
  "High family income", "Medium family income", "Low family income", "Lower educational level",
  "High educational level", "Medium educational level", "Low educational level", "Pre-stroke WHO-5 score",
  "Pre-stroke mRS > 0"
)

df_mmrm_coef_clean <- df_mmrm_coef_clean[!df_mmrm_coef_clean$var %in% c("Lower family income", "Lower educational level"), ]

## Plot
p_list <- lapply(list("all" = list(), "sign" = list(only_sign = TRUE)), \(.x){
  # browser()
  do.call(
    forest_plot,
    modifyList(.x, list(data = df_mmrm_coef_clean))
  )
})


purrr::imap(p_list[2], \(.x, .i){
  ggplot2::ggsave(
    filename = glue::glue("mmrm_coef_{.i}.png"),
    plot = .x,
    height = 14,
    width = 14,
    units = "cm",
    dpi = 300,
    scale = .8
  )
})


############ Prediction model

df_coefs_raw <- get_coefs(path = here::here("1 PA Decline/Fra DDV/241002/mmrm.docx"), index.table = 1) |>
  dplyr::filter(variable != "(Intercept)") |>
  dplyr::select(variable, hop_median, drop_median) |>
  setNames(c("variable", "INCREASE", "DECREASE"))



df_coefs_raw <- get_coefs(path = here::here("1 PA Decline/Fra DDV/240624/pa_change_analyses.docx"), index.table = 2) |>
  dplyr::filter(variable != "(Intercept)") |>
  dplyr::select(variable, hop_median, drop_median) |>
  setNames(c("variable", "INCREASE", "DECREASE"))

## Real work
df_coefs <- df_coefs_raw |>
  dplyr::mutate(dplyr::across(
    tidyselect::all_of(c("INCREASE", "DECREASE")),
    function(.x) {
      # signif(
      as.numeric(.x) # ,
      #   3
      # )
    }
  )) |>
  dplyr::mutate(variable = dplyr::if_else(variable == "Alcohol consumption above recommendations",
    "High alcohol consumption", variable
  )) |>
  ## Important step to keep the data ordered for ggplot
  (function(.y) {
    .y |> dplyr::mutate(variable = factor(variable, levels = rev(.y$variable)))
  })()

## Highest ORs
ls_max <- list(df_coefs[c(1, 2)], df_coefs[c(1, 3)]) |>
  setNames(names(df_coefs)[2:3]) |>
  purrr::map(function(.x) {
    .x |>
      setNames(c("var", "val")) |>
      dplyr::mutate(sorting = abs(log(val))) |>
      dplyr::arrange(1 - sorting) |>
      dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ signif(.x, 2))) |>
      head(5)
  })


p_ls_a <- lapply(ls_max, simple_shap_plot)


purrr::imap(p_ls_a, \(.x, .i){
  ggplot2::ggsave(
    filename = glue::glue("pred_coef_{.i}.png"),
    plot = .x,
    height = 8,
    width = 14,
    units = "cm",
    dpi = 300,
    scale = 1
  )
})



################################################################################
############ Study 2
################################################################################


ls_2_supp <- project.aid::docx2list(path = "/Users/au301842/Library/CloudStorage/OneDrive-Personal/Research/PhD/2 TALOS opfølgning/Manuskript/Supplementary material_v2_0.docx")

pase_0 <- ls_2_supp[["30"]][4:7, -2]

pase_4 <- ls_2_supp[["30"]][10:13, -2]

ls_2_main <- project.aid::docx2list(path = "/Users/au301842/Library/CloudStorage/OneDrive-Personal/Research/PhD/2 TALOS opfølgning/Manuskript/SR/Manuscript_2.docx")

pase_change <- ls_2_main[["169"]][2:5, -2]


pase_clean <- list(
  pase_change,
  pase_0,
  pase_4
) |>
  lapply(pase_table_clean)

pase_clean_plots <- purrr::map2(
  pase_clean,
  list(
    list(
      rev_y = TRUE,
      group.colors = c(
        "#30123BFF",
        "#FABA39FF",
        "#1AE4B6FF",
        "#7A0403FF"
      )
    ),
    list(
      rev_y = FALSE,
      group.colors = viridisLite::viridis(6, option = "B")[2:5]
    ),
    list(
      rev_y = FALSE,
      group.colors = viridisLite::viridis(6, option = "B")[2:5]
    )
  ), \(.x, .y){
    do.call(
      forest_plot_grp,
      modifyList(
        list(
          data = .x
        ),
        .y
      )
    )
  }
)

pase_clean_plots[2]

purrr::imap(setNames(
  pase_clean_plots,
  c(
    "pase_change",
    "pase_0",
    "pase_4"
  )
), \(.x, .i){
  ggplot2::ggsave(
    filename = glue::glue("pase_risk_{.i}.png"),
    plot = .x,
    height = 10,
    width = 12,
    units = "cm",
    dpi = 300,
    scale = 1.2
  )
})

############ Only mRS 0
############ 

pase_change_mrs0_clean <- ls_2_supp[["39"]][3:6,-2] |> pase_table_clean()

pase_change_mrs0_plot <- forest_plot_grp(pase_change_mrs0_clean,
rev_y = TRUE,
group.colors = c(
  "#30123BFF",
  "#FABA39FF",
  "#1AE4B6FF",
  "#7A0403FF"
))

ggplot2::ggsave(
  filename = "pase_risk_mrs0.png",
  plot = pase_change_mrs0_plot,
  height = 10,
  width = 12,
  units = "cm",
  dpi = 300,
  scale = 1.2
)

################################################################################
############ Study 3
################################################################################

ls_3_main <- project.aid::docx2list(path = "/Users/au301842/Library/CloudStorage/OneDrive-Personal/Research/PhD/3 SVD and stroke/Manuskript/EJN/SVD burden_v4_0.docx")
svd <- ls_3_main[["138"]][3:6, ]

svd_clean <- purrr::imap(svd[2:4], \(.x, .i){
  # browser()
  out <- data.frame(or_ci = .x) |>
    tidyr::separate(
      col = "or_ci", into = c("val", "ci"), sep = " \\("
    ) |>
    tidyr::separate(
      col = "ci", into = c("lo", "hi"), sep = " to "
    ) |>
    apply(2, \(.s){
      gsub("\\)", "", .s)
    }) |>
    as.data.frame() |>
    dplyr::mutate(dplyr::across(tidyselect::all_of(c("val", "lo", "hi")), as.numeric))

  out[1, 1] <- 1
  # browser()
  out$var <- sapply(svd[[1]], \(.s){
    l <- nchar(.s)
    substr(.s, 5, l)
  })
  # out[1,2:4] <- 1
  out$model <- .i
  # out$var <- sapply(out$var,\(.x){
  #   paste(strwrap(.x, 25), collapse = "\n")})

  out
}) |> dplyr::bind_rows()

svd_plot <- svd_clean |> 
  forest_plot_grp(rev_y = FALSE, group.colors = viridisLite::viridis(6, option = "B")[2:5])

ggplot2::ggsave(
  filename = "svd_olr_coef.png",
  plot = svd_plot,
  height = 10,
  width = 12,
  units = "cm",
  dpi = 300,
  scale = 1.2
)

# p_colors <- lapply(LETTERS[1:8], \(.x){
#   viridisLite::viridis(6, option = .x)[2:5]
# })
# project.aid::color_plot(p_colors[[2]])
# 
# lapply(p_colors, project.aid::color_plot)