data_preprocessing.Rmd

---
title: "Prepare Harmonized Raw Data for Analysis"
author: "Matteo Delucchi"
output:
  rmarkdown::html_document:
    toc: yes
    toc_float:
      collapsed: true
      smooth_scroll: false
    toc_depth: 3
  rmarkdown::html_vignette: default
vignette: >
  %\VignetteIndexEntry{Prepare Harmonized Raw Data for Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r development, eval=FALSE, include=FALSE}
devtools::document()
devtools::load_all()
renv::snapshot(prompt = F)
```

```{r setup}
# Clear working environment
rm(list=ls())

# Load libraries
library(bnaiaR)
library(tidyr)
library(dplyr)
library(ggplot2)
library(lubridate)
library(RMariaDB)
library(forcats)

# Save plots as files? 
# Warning: Bad layout with knitr. Works well if chunks are run without rendering.
SAVEPLOTS <- FALSE
PLOTPATH <- Sys.getenv("PLOTPATH")
```

# Load Harmonized Raw Data

```{r}
load(system.file("extdata", "adb_raw.RData", package = "bnaiaR"))
str(adb.raw)
```

# Transform Variables: Merging levels, binarize and logarithms

## Remove NAs

```{r}
# NA overview: Filter all, with at least one NA in any column.
plt <- adb.raw %>%
  summarise(across(.cols = everything(), .fns = function(x){is.na(x)})) %>%
  filter(if_any()) %>%
  summarise(across(everything(), sum)) %>%
  
  t() %>%
  as.data.frame() %>%
  mutate(names = rownames(.)) %>%

  ggplot()+
  aes(y=V1, x=names)+
  geom_col() +
  ylab("no. NAs")+
  xlab(element_blank())+
  labs(title = "Missing records in raw data") +
  theme_bw()+
  coord_flip() 

if (SAVEPLOTS) {
  ggsave(
    plot = plt,
    filename = "na_overview_plt.png",
    path = PLOTPATH,
    dpi = 300
  )
} else {
  plt
}
```


```{r}
length(adb.raw$Age.at.diag.grouped[which(adb.raw$Age.at.diag.grouped == "H" | is.na(adb.raw$Age.at.diag.grouped))])
```


```{r}
df <- adb.raw %>% 
  na.omit() 
str(df)
```


```{r, fig.height=14, out.height= "100%"}
plt <- overview_plt_max(df)

if (SAVEPLOTS){
  ggsave(plot = plt, filename = "ADB_noNAs_overview_raw.png", path = PLOTPATH, 
         width = 14, height = 20, dpi = 300)
  }
```

## Rename Variables

```{r}
df <- df %>% 
  rename(Hypertension = Hypertension) %>%
  rename(AgeDiag = age.at.time.of.diagnosis) %>%
  rename(AgeDiag.group = Age.at.diag.grouped) %>%
  rename(IAsize = IA_size.at.time.of.SAH) %>%
  rename(IAsize.groups = IA_size.at.time.of.SAH.groups) %>%
  rename(Basis.of.recruitement = Basis.of.recruitment_CaseControl)

str(df)
```

## Large Age groups reduce number of levels

```{r}
plt.age <- df %>% 
  # add layouting information
  group_by(AgeDiag) %>%
  summarise(n = n()) %>%
  mutate(label = case_when(n < 5 ~ n)) %>%
  
  # Plot
  ggplot()+
  aes(x=AgeDiag, y=n)+
  geom_col(width = 0.5) +
  geom_text(aes(label = label, vjust = -0.2))+
  theme_minimal()+
  labs(title= "Age at time of Diagnosis",
          # caption = stringr::str_wrap("Age at time of diagnosis of IA."),
       y="Number of records",
       x="Age at time of IA diagnosis")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")

plt.age
```

Age at time of diagnosis is discretized following clinical conventions and
grouped in bins of 5y between the age of 40 and 65.

```{r}
plt.age5y <- df %>% 
  # add layouting information
  group_by(AgeDiag.group) %>%
  summarise(n = n()) %>%
  mutate(label = case_when(n < 250 ~ n)) %>%
  
  # Plot
  ggplot()+
  aes(x=AgeDiag.group, y=n)+
  geom_col() +
  geom_text(aes(label = label, vjust = -0.2))+
  theme_minimal()+
  labs(title= "Case-distribution of grouped Age at time of Diagnosis",
          subtitle = "Grouped in bins of 5y for patients between 40y and 65y.",
          caption = stringr::str_wrap("Age at time of diagnosis of IA grouped in A=[0,40), B=[40,45), C=[45,50), D=[50,55), E=[55,60), F=[60,65), G=[65, max. AgeDiag)"),
       y="Number of records",
       x="Age Group")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")

if (SAVEPLOTS){
  ggsave(plot = plt.age5y, filename = "age_grouped_5y.png", path = PLOTPATH, dpi = 600)
} else {
  plt.age5y
}
```


To reduce the number of classes of the AgeDiag.grouped variable a coarser 
grouping is included.

```{r}
df <- df %>%
    mutate(
      AgeDiag.group.coarse = cut(
        AgeDiag,
        breaks = c(0,
                   seq(40, 60, 10),
                   max(AgeDiag, na.rm = TRUE)),
        include.lowest = TRUE,
        right = FALSE
      )
    ) %>%
    mutate(AgeDiag.group.coarse = factor(AgeDiag.group.coarse, labels = LETTERS[1:length(levels(AgeDiag.group.coarse))])) 
```


```{r}
plt.age10y <- df %>%
  # add layouting information
  group_by(AgeDiag.group.coarse) %>%
  summarise(n = n()) %>%
  mutate(label = case_when(n < 300 ~ n)) %>%
  
  # Plot
  ggplot()+
  aes(x=AgeDiag.group.coarse, y=n)+
  geom_col() +
  geom_text(aes(label = label, vjust = -0.2))+
  theme_minimal()+
  labs(title= "Case-distribution of grouped Age at time of Diagnosis",
          subtitle = "Grouped in bins of 10y for patients between 40y and 60y.",
          caption = stringr::str_wrap("Age at time of diagnosis of IA grouped in A=[0,40), B=[40,50), C=[50,60), D=[60,max. AgeDiag)"),
       y="Number of records",
       x="Age Group")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")

if (SAVEPLOTS){
  ggsave(plot = plt.age10y, filename = "age_grouped_10y.png", path = PLOTPATH, dpi = 600)
} else {
  plt.age10y
}
```

Multiplot of age

```{r, out.height="100%", out.width="100%"}
plt.age.mult <- plt.age+
  labs(title= "",
          subtitle = "",
          caption = "",
       y="Number of records",
       x="Age")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA))

plt.age5y.mult <- plt.age5y +
  labs(title= "",
          subtitle = "",
          caption = "",
       y="Number of records",
       x="Age Group")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA))

plt.age10y.mult <- plt.age10y+
  labs(title= "",
          subtitle = "",
          caption = "",
       y="Number of records",
       x="Age Group")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA))

plt.comb.age <- cowplot::plot_grid(plt.age.mult, plt.age5y.mult, plt.age10y.mult, labels = "AUTO", ncol = 1)+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA))

if (SAVEPLOTS){
  ggsave(plot = plt.comb.age, filename = "age_multiplot.png", path = PLOTPATH, 
         dpi = 600, height = 13, bg = "transparent")
} else {
  plt.comb.age
}
```
Overlay plot age

```{r}
df.age <- df %>% 
  select(c(AgeDiag, AgeDiag.group, AgeDiag.group.coarse)) %>%
  arrange(AgeDiag)

df.age.cont <- df %>% 
  select(c(AgeDiag, AgeDiag.group, AgeDiag.group.coarse)) %>%
  arrange(AgeDiag) %>%
  count(AgeDiag, name="AgeDiag.n")

df.age.group <- df %>% 
  select(c(AgeDiag, AgeDiag.group, AgeDiag.group.coarse)) %>%
  arrange(AgeDiag) %>%
  count(AgeDiag.group,  name="AgeDiag.group.n")

df.age.group.coarse <- df %>% 
  select(c(AgeDiag, AgeDiag.group, AgeDiag.group.coarse)) %>%
  arrange(AgeDiag) %>%
  count(AgeDiag.group.coarse, name="AgeDiag.group.coarse.n")
```


```{r}
plt.overlay.age <- df.age.cont %>%
  ggplot() +
  # Age group coarse
  annotate(
    "rect",
    xmin = c(
      min(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "A")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "B")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "C")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "D")])
    ),
    xmax = c(
      max(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "A")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "B")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "C")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group.coarse == "D")])
    ),
    ymin = rep(0, 4),
    ymax = df.age.group.coarse$AgeDiag.group.coarse.n,
    alpha = 0.2,
    color = "orange",
    fill = "yellow"
  ) +
  # Age group fine
  annotate(
    "rect",
    xmin = c(
      min(df.age$AgeDiag[which(df.age$AgeDiag.group == "A")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group == "B")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group == "C")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group == "D")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group == "E")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group == "F")]),
      min(df.age$AgeDiag[which(df.age$AgeDiag.group == "G")])
    ),
    xmax = c(
      max(df.age$AgeDiag[which(df.age$AgeDiag.group == "A")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group == "B")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group == "C")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group == "D")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group == "E")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group == "F")]),
      max(df.age$AgeDiag[which(df.age$AgeDiag.group == "G")])
    ),
    ymin = rep(0, 7),
    ymax = df.age.group$AgeDiag.group.n,
    alpha = 0.2,
    color = "blue",
    fill = "blue"
  ) +
  # Age group cont
  geom_col(aes(x = AgeDiag, y = AgeDiag.n), width = 0.5) +
  # Layouting
  theme_minimal() +
  scale_x_continuous(breaks = c(0, 40, seq(40, 65, 10), max(df.age.cont$AgeDiag)),
                     minor_breaks = seq(40, 65, 5)) +
  labs(
    title = "Case-distribution of age at time of IA diagnosis",
    subtitle = "Grouped in bins of 5y (blue) and 10y (yellow) for patients between 40y and 65y.",
    y = "Number of records",
    x = "Age at diagnosis"
  ) +
  theme(
    plot.caption = element_text(hjust = 0.5),
    plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA)
  )

if (SAVEPLOTS) {
  ggsave(
    plot = plt.overlay.age,
    filename = "age_overlay.png",
    path = PLOTPATH,
    dpi = 600
  )
} else {
  plt.overlay.age
}
```


## Binarize Age groups

```{r}
AgeDiag.group_bin <- df %>% 
  select(AgeDiag.group) %>% 
  correlationfunnel::binarize(thresh_infreq = 0.0001)

df <- df %>%
  cbind(AgeDiag.group_bin) %>%
  mutate(across(contains("AgeDiag.group__"), as.factor))
  
str(df)
```



## Merge and Binarize Smoking

```{r}
Smoking_Current_Former_No_bin <- df %>% 
  select(Smoking_Current_Former_No) %>% 
  correlationfunnel::binarize(thresh_infreq = 0.0001)

df <- df %>%
  cbind(Smoking_Current_Former_No_bin) %>%
  mutate(across(contains("Smoking_Current_Former_No__"), as.factor))
  
str(df)
```

Merging former and non-smokers in "NotCurrent" smokers:

```{r}
df <- df %>%
  # Make smoking binary:
  mutate(Smoking_Current_Former_No = as.character(Smoking_Current_Former_No)) %>%
  mutate(Smoking_CurrentNotCurrent = case_when(Smoking_Current_Former_No == "Current" ~ "Current",
                                               Smoking_Current_Former_No == "Former" ~ "NotCurrent",
                                               Smoking_Current_Former_No == "No" ~ "NotCurrent",
                                               TRUE ~ Smoking_Current_Former_No)) %>% 
    mutate(Smoking_Current_Former_No = factor(Smoking_Current_Former_No),
           Smoking_CurrentNotCurrent = factor(Smoking_CurrentNotCurrent)) 

str(df)
```


## Creation of the cohort of patients with lesions in different risk locations

We group IA location by vessel into three classes according their risk of rupture:

```{r}
df <- df %>%
  mutate(location.grouped = if_else(IA_Location == "Acom" |
                          IA_Location == "Pcom" |
                          IA_Location == "V-B" |
                          IA_Location == "A2" |
                          IA_Location == "PC", 
                                       "High", 
                                       if_else(IA_Location == "MCA" |
                                                 IA_Location == "ICA" |
                                                 IA_Location == "Basilar" |
                                                 IA_Location == "Other" |
                                                 IA_Location == "A1 segment ant",
                                               "Medium",
                                               if_else(IA_Location == "CavICA" |
                                                         IA_Location == "OphtICA",
                                                       "Low",
                                                       NA_character_))))%>%
  mutate_at(vars(location.grouped), list(~as.factor(.)))
str(df)
```

And additionally we binarize all locations

```{r}
IA_Location_bin <- df %>% 
  select(IA_Location) %>% 
  correlationfunnel::binarize(thresh_infreq = 0.0001)%>%
  mutate(IA_Location__VB = `IA_Location__V-B`) %>%
  select(-c(`IA_Location__V-B`)) 
IA_Location_bin <- IA_Location_bin %>%
  select(sort(names(IA_Location_bin)))

location.grouped <- df %>% 
  select(location.grouped) %>% 
  correlationfunnel::binarize(thresh_infreq = 0.0001)

df <- df %>%
  cbind(IA_Location_bin) %>%
  cbind(location.grouped) %>%
  mutate(across(contains("IA_Location__"), as.factor),
         across(contains("location.grouped__"), as.factor))
str(df)
```

### Plots

```{r}
df.loc <- df %>%
  select(c(IA_Location, location.grouped))
df.loc
```



```{r}
plt.locv <- df %>%
  # add layouting information
  group_by(IA_Location) %>%
  summarise(n = n()) %>%
  mutate(label = case_when(n < 300 ~ n)) %>%
  
  # Plot
  ggplot()+
  aes(x=IA_Location, y=n)+
  geom_col() +
  geom_text(aes(label = label, hjust = -0.2))+
  theme_minimal()+
  labs(title= "Case-distribution of IA location by vessel",
       y="Number of records",
       x=NULL)+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")+
    coord_flip()

plt.locv
```

```{r}
plt.locr <- df %>%
  # Plot
  ggplot()+
  aes(x=location.grouped)+
  geom_bar() +
  theme_minimal()+
  labs(title= "Case-distribution of IA location by risk of rupture",
       y="Number of records",
       x=NULL)+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")+
    coord_flip()

plt.locr
```

```{r}
plt.loc.comb <- df.loc %>%
  # Plot
  ggplot()+
  aes(x=IA_Location, fill=location.grouped)+
  geom_bar() +
  theme_minimal()+
  labs(title= "Case-distribution of IA location",
       y="Number of records",
       x=NULL)+
  scale_fill_manual(values=c("High" = "darkred", "Medium"= "darkorange", "Low"="darkgreen"), name = "Risk of IA rupture")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")+
  coord_flip()

if (SAVEPLOTS) {
  ggsave(
    plot = plt.loc.comb,
    filename = "location_multiplot.png",
    path = PLOTPATH,
    dpi = 600
  )
} else {
  plt.loc.comb
}
```

## Size groups merged to reduce number of levels

The two classes of large IAs have a small number of observations.
We merge them.

```{r}
df <- df %>%
  # merge size groups to reduce number of levels
  mutate(IAsize.groups = as.character(IAsize.groups)) %>%
  mutate(IAsize.groups.merged = case_when(IAsize.groups == "D" ~ "C",
                                                          TRUE ~ IAsize.groups)) %>%
  filter(IAsize.groups != "E")%>%
  mutate(IAsize.groups = as.factor(IAsize.groups),
         IAsize.groups.merged = as.factor(IAsize.groups.merged))
```

And additionally we binarize all size groups

```{r}
IA_size_groups_bin <- df %>% 
  select(IAsize.groups) %>% 
  correlationfunnel::binarize(thresh_infreq = 0.0001)

IA_size_groups_merged_bin <- df %>% 
  select(IAsize.groups.merged) %>% 
  correlationfunnel::binarize(thresh_infreq = 0.0001)

df <- df %>%
  cbind(IA_size_groups_bin) %>%
  cbind(IA_size_groups_merged_bin) %>%
  mutate(across(contains("IAsize.groups__"), as.factor)) %>%
  mutate(across(contains("IAsize.groups.merged__"), as.factor))

str(df)
```

```{r}
plt.size.comb <- df %>%
  # Plot
  ggplot()+
  aes(x=IAsize.groups, fill=IAsize.groups.merged)+
  geom_bar() +
  theme_minimal()+
  labs(title= "Case-distribution of IA size groups",
       y="Number of records",
       x=NULL)+
  scale_fill_manual(values=c("A"=palette.colors(4)[[2]], "B"=palette.colors(4)[[3]], "C"=palette.colors(4)[[4]]),
                    name = "IA size groups merged",
                    labels = LETTERS[1:3])+
  coord_flip()+
  theme(legend.position = c(0.8,0.8),
        legend.background = element_rect(fill = "white", color= "gray"))
plt.size.comb
```

## Size logarithmic fits gaussian distribution


```{r}
df <- df %>%
  # add logarithm of size
  mutate(IAsize_log = log(round(IAsize, 2))) %>%
  filter(IAsize_log > 0)

str(df)
```


IAsize outliers

```{r}
plt.size.log <- ggplot(df, aes(x= IAsize_log)) +
      geom_bar(stat = "count", show.legend = TRUE, col="black", fill="gray") +
      geom_text(aes(label=ifelse((stat(count)<10) & (x>4), stat(count),"")), stat='count', nudge_y=20)+
      theme_minimal()+
  labs(title= "Case-distribution of log(IAsize)",
          caption = stringr::str_wrap("Histogram of log(IAsize) with large values (log(IAsize)>4) labeled with their number of records."),
       y="Number of records",
       x="log(IAsize)")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")

if (SAVEPLOTS){
  ggsave(plot = plt.size.log, filename = "IAsize_log_outliers.png", path = PLOTPATH,
         width =9, height = 7, dpi = 600)
} else {
  plt.size.log
}
```

Multiplot IA size

```{r}
# IA size cont.
plt.size <- ggplot(df, aes(x= IAsize)) +
      geom_bar(stat = "count", show.legend = TRUE, col="black", fill="gray") +
      # geom_text(aes(label=ifelse((stat(count)<10) & (x>4), stat(count),"")), stat='count', nudge_y=20)+
      theme_minimal()+
  scale_x_continuous(breaks = c(seq(0, 80, 20), seq(80, max(df$IAsize), 50)),
                     minor_breaks = seq(80, max(df$IAsize), 50))+
  labs(title= "Case-distribution of IA size",
          # caption = stringr::str_wrap("Histogram of IA size with large values (IAsize>50) labeled with their number of records."),
       y="Number of records",
       x="max. Diameter [mm]")+
  theme(plot.caption = element_text(hjust = 0.5),
        plot.caption.position = "plot")

# IA size log.
plt.size.log <- plt.size.log+
  theme(plot.caption = element_blank(),
        plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA))

# IA size grouped
plt.size.comb <- plt.size.comb +
  theme(plot.caption = element_blank(),
        plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA))

plt.comb.size <- cowplot::plot_grid(plt.size, plt.size.log, plt.size.comb, labels = "AUTO", ncol = 1)+
  theme(plot.caption = element_blank(),
        plot.caption.position = "plot",
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", colour = NA))

if (SAVEPLOTS){
  ggsave(plot = plt.comb.size, filename = "size_multiplot.png", path = PLOTPATH, 
         dpi = 600, height = 13, bg = "transparent")
} else {
  plt.comb.size
}

```


## Order levels in "normal" order
```{r}
df <- df %>%
  mutate(Positive.famillial.history = fct_relevel(Positive.famillial.history, c("Yes","No")))%>%
  mutate(Smoking_Current_Former_No = fct_relevel(Smoking_Current_Former_No, c("Current","Former","No")))%>%
  mutate(Ruptured_IA = fct_relevel(Ruptured_IA, c("Yes","No")))%>%
  mutate(Multiple.IAs = fct_relevel(Multiple.IAs, c("Yes", "No")))%>%
  mutate(IAsize.groups = fct_relevel(IAsize.groups, LETTERS[1:length(levels(IAsize.groups))]))%>%
  mutate(AgeDiag.group = fct_relevel(AgeDiag.group, LETTERS[1:length(levels(AgeDiag.group))])) %>%
  mutate(Hypertension = fct_relevel(Hypertension, c("AnyType", "Never")))%>%
  mutate(IA_Location = fct_relevel(IA_Location, c("Acom", "Pcom","V-B", "A2","PC","MCA","ICA","Basilar","Other","A1 segment ant","CavICA","OphtICA")))%>%
  mutate(Gender = fct_relevel(Gender, c("Female", "Male")))%>%
  mutate(location.grouped = fct_relevel(location.grouped, c("High", "Medium", "Low")))
str(df)
```

## Save ADB as .rds 

```{r save adb}
adb <- df
str(as.data.frame(adb))

# saveRDS(adb, 
#         file = "/home/delt/ZHAW/BNstructureLearning/2108_adb-v0.2.rds")
```

# Include Preprocessed Data in bnaiaR package

```{r}
# usethis::use_data(adb, overwrite = TRUE)
# devtools::document()
```