AEDB.EM.baseline.Rmd

---
title: "Baseline characteristics"
author: "[Sander W. van der Laan, PhD](https://swvanderlaan.github.io) | @swvanderlaan | s.w.vanderlaan@gmail.com"
date: "`r Sys.Date()`"
output:
  html_notebook:
    cache: yes
    code_folding: hide
    collapse: yes
    df_print: paged
    fig.align: center
    fig_caption: yes
    fig_height: 6
    fig_retina: 2
    fig_width: 7
    highlight: tango
    theme: lumen
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: yes
mainfont: Arial
subtitle: Accompanying 'EntropyMasker'
editor_options:
  chunk_output_type: inline
# bibliography: references.bib
# knit: worcs::cite_all
---

# General Setup
We will clean the environment, setup the locations, define colors, and create a datestamp.

_Clean the environment._
```{r echo = FALSE}
rm(list = ls())
```

_Set locations and working directories..._
```{r LocalSystem, echo = FALSE}
source("scripts/local.system.R")

```

_... a package-installation function ..._
```{r Function: installations}
source("scripts/functions.R")

```


_... and load those packages._
```{r loading_packages, message=FALSE, warning=FALSE}
source("scripts/pack01.packages.R")

```

_We will create a datestamp and define the Utrecht Science Park Colour Scheme_.
```{r Setting: Colors}

Today = format(as.Date(as.POSIXlt(Sys.time())), "%Y%m%d")
Today.Report = format(as.Date(as.POSIXlt(Sys.time())), "%A, %B %d, %Y")

source("scripts/colors.R")

```

```{r global_options, include = FALSE}
# further define some knitr-options.
knitr::opts_chunk$set(fig.width = 12, fig.height = 8, fig.path = 'Figures/', 
                      wwarning = TRUE, # show warnings during codebook generation
  message = TRUE, # show messages during codebook generation
  error = TRUE, # do not interrupt codebook generation in case of errors,
                # usually better for debugging
  echo = TRUE,  # show R code
                      eval = TRUE)
ggplot2::theme_set(ggplot2::theme_minimal())
pander::panderOptions("table.split.table", Inf)
```

# This notebook 

In this notebook we create a baseline table of the samples used in **`EntropyMasker`**. 

# Athero-Express Biobank Study

The [*Athero-Express Biobank Study (AE)*](https://doi.org/10.1007/s10564-004-2304-6) contains plaque material of patients that underwent endarterectomyat two Dutch tertiary referral centers. Details of the study design were described before. Briefly, blood and plaque material were obtained during endarterectomy and stored at -80 ℃. All patients provided informed consent and the study was approved by the medical ethics committee.

## Load relevant samples

```{r}
EM_samples <- fread(paste0(ANALYSIS_loc, "/dataverse/EntropyMasker_image_files_used.txt"))
```

## Load data

Loading Athero-Express Biobank Study clinical and biobank data.

```{r LoadAEDB}
cat("* get Athero-Express Biobank Study Database...")
# METHOD 1: It seems this method gives loads of errors and warnings, which all are hard to comprehend
#           or debug. We expect 3,527 samples, and 927 variables; we get 927 variables!!!
# AEdata = as.data.table(read.spss(paste0(INP_loc,"/2017-1NEW_AtheroExpressDatabase_ScientificAE_20171306_v1.0.sav"),
#                                  trim.factor.names = TRUE, trim_values = TRUE, # we trim spaces in values
#                                  reencode = TRUE, # we re-encode to the local locale encoding
#                                  add.undeclared.levels = "append", # we do *not* want to convert to R-factors
#                                  use.value.labels = FALSE, # we do *not* convert variables with value labels into R factors
#                                  use.missings = TRUE, sub = "NA", # we will set every missing variable to NA
#                                  duplicated.value.labels = "condense", # we will condense duplicated value labels
#                                  to.data.frame = TRUE))
# AEdata.labels <- as.data.table(attr(AEdata, "variable.labels"))
# names(AEdata.labels) <- "Variable"

# METHOD 2: Using library("haven") importing seems flawless; best argument being:
#           we expect 3,527 samples and 888 variables, which is what you'd get with this method
#           So for now, METHOD 2 is prefered. 
#            
require(haven)

# AEDB <- haven::read_sav(paste0(AEDB_loc, "/2022_1_NEW_AtheroExpressDatabase_ScientificAE_15-02-2022.sav")) # something wrong with Age-variable
# AEDB <- haven::read_sav(paste0(AEDB_loc, "/2020_1_NEW_AtheroExpressDatabase_ScientificAE_30-09-2020.sav")) # duplicate studynumbers in it
AEDB <- haven::read_sav(paste0(AEDB_loc, "/2020_1_NEW_AtheroExpressDatabase_ScientificAE_16-03-2020.sav"))

# writing off the SPSS data to an Excel.
# fwrite(AEdata, file = paste0(INP_loc,"/2017-1NEW_AtheroExpressDatabase_ScientificAE_20171306_v1.0.values.xlsx"), 
#        sep = ";", na = "NA", dec = ".", col.names = TRUE, row.names = FALSE,
#        dateTimeAs = "ISO", showProgress = TRUE, verbose = TRUE)
# warnings()

AEDB[1:10, 1:10]
dim(AEDB)

```



<!-- ## Examine AEDB -->

<!-- We can examine the contents of the Athero-Express Biobank dataset to know what each variable is called, what class -->
<!-- (type) it has, and what the variable description is. -->

<!-- There is an excellent post on this: <https://www.r-bloggers.com/working-with-spss-labels-in-r/>. -->

<!-- ```{r AEDB: describe} -->
<!-- AEDB %>% sjPlot::view_df(show.type = TRUE, -->
<!--                          show.frq = TRUE, -->
<!--                          show.prc = TRUE, -->
<!--                          show.na = TRUE,  -->
<!--                          max.len = TRUE,  -->
<!--                          wrap.labels = 20, -->
<!--                          verbose = FALSE,  -->
<!--                          use.viewer = FALSE, -->
<!--                          file = paste0(OUT_loc, "/", Today, ".AEDB.dictionary.html"))  -->
<!-- ``` -->

## Fixing and creating variables

We have to fix certain clinical parameters:

-   symptoms
-   diabetes
-   alcohol use
-   smoking
-   plaque phenotypes

### Symptoms

We need to be very strict in defining *symptoms.* Therefore we will fix a new
variable that groups *symptoms* at inclusion.

Coding of *symptoms* is as follows:

-   missing -999
-   Asymptomatic 0
-   TIA 1
-   minor stroke 2
-   Major stroke 3
-   Amaurosis fugax 4
-   Four vessel disease 5
-   Vertebrobasilary TIA 7
-   Retinal infarction 8
-   Symptomatic, but aspecific symtoms 9
-   Contralateral symptomatic occlusion 10
-   retinal infarction 11
-   armclaudication due to occlusion subclavian artery, CEA needed for bypass
    12
-   retinal infarction + TIAs 13
-   Ocular ischemic syndrome 14
-   ischemisch glaucoom 15
-   subclavian steal syndrome 16
-   TGA 17

We will group as follows:

1.  Asymptomatic > 0
2.  TIA > 1, 7, 13
3.  Stroke > 2, 3
4.  Ocular > 4, 14, 15
5.  Retinal infarction > 8, 11
6.  Other > 5, 9, 10, 12, 16, 17

```{r FixSymptoms, message=FALSE, warning=FALSE}

# Fix symptoms

attach(AEDB)
AEDB[,"Symptoms.5G"] <- NA
AEDB$Symptoms.5G[sympt == 0] <- "Asymptomatic"
AEDB$Symptoms.5G[sympt == 1 | sympt == 7 | sympt == 13] <- "TIA"
AEDB$Symptoms.5G[sympt == 2 | sympt == 3] <- "Stroke"
AEDB$Symptoms.5G[sympt == 4 | sympt == 14 | sympt == 15 ] <- "Ocular"
AEDB$Symptoms.5G[sympt == 8 | sympt == 11] <- "Retinal infarction"
AEDB$Symptoms.5G[sympt == 5 | sympt == 9 | sympt == 10 | sympt == 12 | sympt == 16 | sympt == 17] <- "Other"


# AsymptSympt
AEDB[,"AsymptSympt"] <- NA
AEDB$AsymptSympt[sympt == -999] <- NA
AEDB$AsymptSympt[sympt == 0] <- "Asymptomatic"
AEDB$AsymptSympt[sympt == 1 | sympt == 7 | sympt == 13 | sympt == 2 | sympt == 3] <- "Symptomatic"
AEDB$AsymptSympt[sympt == 4 | sympt == 14 | sympt == 15 | sympt == 8 | sympt == 11 | sympt == 5 | sympt == 9 | sympt == 10 | sympt == 12 | sympt == 16 | sympt == 17] <- "Ocular and others"

# AsymptSympt
AEDB[,"AsymptSympt2G"] <- NA
AEDB$AsymptSympt2G[sympt == -999] <- NA
AEDB$AsymptSympt2G[sympt == 0] <- "Asymptomatic"
AEDB$AsymptSympt2G[sympt == 1 | sympt == 7 | sympt == 13 | sympt == 2 | sympt == 3 | sympt == 4 | sympt == 14 | sympt == 15 | sympt == 8 | sympt == 11 | sympt == 5 | sympt == 9 | sympt == 10 | sympt == 12 | sympt == 16 | sympt == 17] <- "Symptomatic"

detach(AEDB)

# table(AEDB$sympt, useNA = "ifany")
# table(AEDB$AsymptSympt2G, useNA = "ifany")
# table(AEDB$Symptoms.5G, useNA = "ifany")
# 
# table(AEDB$AsymptSympt2G, AEDB$sympt, useNA = "ifany")
# table(AEDB$Symptoms.5G, AEDB$sympt, useNA = "ifany")
table(AEDB$AsymptSympt2G, AEDB$Symptoms.5G, useNA = "ifany")

# AEDB.temp <- subset(AEDB,  select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "sympt", "Symptoms.5G", "AsymptSympt"))
# require(labelled)
# AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
# AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
# AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
# 
# DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
# 
# table(AEDB.temp$Symptoms.5G, AEDB.temp$AsymptSympt)
# 
# rm(AEDB.temp)
```

### Re-assessed symptoms

We re-assessed the categorization of symptoms. These are summarized and parsed
in this section.

Labeling of new symptom categories.

```{r message=FALSE, warning=FALSE, paged.print=TRUE}

AEDB$indexsymptoms_worst
AEDB$indexsymptoms_worst_4g
AEDB$indexsymptoms_latest
AEDB$indexsymptoms_latest_4g

```

Getting counts for each of the most important categories.

```{r}
cat("New 'worst' vs 'latest' symptom categories.\n")
table(as_factor(AEDB$indexsymptoms_worst_4g), as_factor(AEDB$indexsymptoms_latest_4g))

cat("\nNew 'worst' symptom categories.\n")
table((AEDB$indexsymptoms_worst_4g))

cat("\nNew 'latest' symptom categories.\n")
table(as_factor(AEDB$indexsymptoms_latest_4g))
```

Comparing with the original symptom categories.

```{r}
cat("New 'latest' vs original symptom 2G categories.\n")
table((AEDB$indexsymptoms_latest_4g), AEDB$AsymptSympt2G)

cat("\nNew 'latest' vs original symptom 5G categories.\n")
table((AEDB$indexsymptoms_latest_4g), AEDB$Symptoms.5G)
   
```

We need to be very strict in defining *symptoms.* Therefore we will fix a new
variable that groups *symptoms* at inclusion.

Coding of *symptoms* is as follows:

-   asympt 0\
-   ocular 1\
-   TIA 2\
-   stroke 3\
-   unclear 9

We will group as follows:

1.  Asymptomatic > 0
2.  Symptomatic > 1, 2, 3
3.  NA > 9

```{r FixNewSymptoms}

# Fix symptoms
attach(AEDB)

# Symptoms.Update2G
AEDB[,"Symptoms.Update2G"] <- NA
AEDB$Symptoms.Update2G[indexsymptoms_latest_4g == 0] <- "Asymptomatic"
AEDB$Symptoms.Update2G[indexsymptoms_latest_4g == 1 | indexsymptoms_latest_4g == 2 | indexsymptoms_latest_4g == 3] <- "Symptomatic"
AEDB$Symptoms.Update2G[indexsymptoms_latest_4g == 9 ] <- NA

# Symptoms.Update3G
AEDB[,"Symptoms.Update3G"] <- NA
AEDB$Symptoms.Update3G[indexsymptoms_latest_4g == 0] <- "Asymptomatic"
AEDB$Symptoms.Update3G[indexsymptoms_latest_4g == 1 | indexsymptoms_latest_4g == 2 | indexsymptoms_latest_4g == 3] <- "Symptomatic"
AEDB$Symptoms.Update3G[indexsymptoms_latest_4g == 9 ] <- "Unclear"

detach(AEDB)

table(AEDB$Symptoms.Update2G, AEDB$Symptoms.5G, useNA = "ifany")
table(AEDB$Symptoms.Update3G, AEDB$Symptoms.5G, useNA = "ifany")

```

### Other clinical characteristics

We will also fix the *diabetes* status variable.

```{r FixDiabetes, message=FALSE, warning=FALSE}

# Fix diabetes
attach(AEDB)
AEDB[,"DiabetesStatus"] <- NA
AEDB$DiabetesStatus[DM.composite == -999] <- NA
AEDB$DiabetesStatus[DM.composite == 0] <- "Control (no Diabetes Dx/Med)"
AEDB$DiabetesStatus[DM.composite == 1] <- "Diabetes"
detach(AEDB)

table(AEDB$DM.composite, AEDB$DiabetesStatus)
# AEDB.temp <- subset(AEDB,  select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "DM.composite", "DiabetesStatus"))
# require(labelled)
# AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
# AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
# AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
# AEDB.temp$DiabetesStatus <- to_factor(AEDB.temp$DiabetesStatus)
# 
# DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
# 
# rm(AEDB.temp)

```

We will also fix the *smoking* status variable. We are interested in whether
someone never, ever or is currently (at the time of inclusion) smoking. This is
based on the questionnaire.

-   `diet801`: are you a smoker?
-   `diet802`: did you smoke in the past?

We already have some variables indicating smoking status:

-   `SmokingReported`: patient has reported to smoke.
-   `SmokingYearOR`: smoking in the year of surgery?
-   `SmokerCurrent`: currently smoking?

```{r FixSmoking, message=FALSE, warning=FALSE}
require(labelled)
AEDB$diet801 <- to_factor(AEDB$diet801)
AEDB$diet802 <- to_factor(AEDB$diet802)
AEDB$diet805 <- to_factor(AEDB$diet805)
AEDB$SmokingReported <- to_factor(AEDB$SmokingReported)
AEDB$SmokerCurrent <- to_factor(AEDB$SmokerCurrent)
AEDB$SmokingYearOR <- to_factor(AEDB$SmokingYearOR)

# table(AEDB$diet801)
# table(AEDB$diet802)
# table(AEDB$SmokingReported)
# table(AEDB$SmokerCurrent)
# table(AEDB$SmokingYearOR)
# table(AEDB$SmokingReported, AEDB$SmokerCurrent, useNA = "ifany", dnn = c("Reported smoking", "Current smoker"))
# 
# table(AEDB$diet801, AEDB$diet802, useNA = "ifany", dnn = c("Smoker", "Past smoker"))

cat("\nFixing smoking status.\n")
attach(AEDB)
AEDB[,"SmokerStatus"] <- NA
AEDB$SmokerStatus[diet802 == "don't know"] <- "Never smoked"
AEDB$SmokerStatus[diet802 == "I still smoke"] <- "Current smoker"
AEDB$SmokerStatus[SmokerCurrent == "no" & diet802 == "no"] <- "Never smoked"
AEDB$SmokerStatus[SmokerCurrent == "no" & diet802 == "yes"] <- "Ex-smoker"
AEDB$SmokerStatus[SmokerCurrent == "yes"] <- "Current smoker"
AEDB$SmokerStatus[SmokerCurrent == "no data available/missing"] <- NA
# AEDB$SmokerStatus[is.na(SmokerCurrent)] <- "Never smoked"
detach(AEDB)

cat("\n* Current smoking status.\n")
table(AEDB$SmokerCurrent,
      useNA = "ifany", 
      dnn = c("Current smoker"))

cat("\n* Updated smoking status.\n")
table(AEDB$SmokerStatus,
      useNA = "ifany", 
      dnn = c("Updated smoking status"))

cat("\n* Comparing to 'SmokerCurrent'.\n")
table(AEDB$SmokerStatus, AEDB$SmokerCurrent, 
      useNA = "ifany", 
      dnn = c("Updated smoking status", "Current smoker"))

# AEDB.temp <- subset(AEDB,  select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "DM.composite", "DiabetesStatus"))
# require(labelled)
# AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
# AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
# AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
# AEDB.temp$DiabetesStatus <- to_factor(AEDB.temp$DiabetesStatus)
# 
# DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
# 
# rm(AEDB.temp)


```

We will also fix the *alcohol* status variable.

```{r FixAlcohol, message=FALSE, warning=FALSE}

# Fix diabetes
attach(AEDB)
AEDB[,"AlcoholUse"] <- NA
AEDB$AlcoholUse[diet810 == -999] <- NA
AEDB$AlcoholUse[diet810 == 0] <- "No"
AEDB$AlcoholUse[diet810 == 1] <- "Yes"
detach(AEDB)

# AEDB.temp <- subset(AEDB,  select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "diet810", "AlcoholUse"))
# require(labelled)
# AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
# AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
# AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
# AEDB.temp$AlcoholUse <- to_factor(AEDB.temp$AlcoholUse)
# 
# DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
# 
# rm(AEDB.temp)


```

We will also fix a history of CAD, stroke or peripheral intervention status variable. This will be based on `CAD_history`, `Stroke_history`, and `Peripheral.interv`

```{r FixCAD_History, message=FALSE, warning=FALSE}

# Fix diabetes
attach(AEDB)
AEDB[,"MedHx_CVD"] <- NA
AEDB$MedHx_CVD[CAD_history == 0 | Stroke_history == 0 | Peripheral.interv == 0] <- "No"
AEDB$MedHx_CVD[CAD_history == 1 | Stroke_history == 1 | Peripheral.interv == 1] <- "yes"
detach(AEDB)

table(AEDB$CAD_history)
table(AEDB$Stroke_history)
table(AEDB$Peripheral.interv)
table(AEDB$MedHx_CVD)

# AEDB.temp <- subset(AEDB,  select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "diet810", "AlcoholUse"))
# require(labelled)
# AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
# AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
# AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
# AEDB.temp$AlcoholUse <- to_factor(AEDB.temp$AlcoholUse)
# 
# DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
# 
# rm(AEDB.temp)


```

### Plaque phenotypes

We will also fix the *plaquephenotypes* variable.

Coding of symptoms is as follows:

-   missing -999\
-   not relevant -888
-   fibrous 1\
-   fibroatheromatous 2\
-   atheromatous 3

```{r FixPlaquePhenotypes, message=FALSE, warning=FALSE}

# Fix plaquephenotypes
attach(AEDB)
AEDB[,"OverallPlaquePhenotype"] <- NA
AEDB$OverallPlaquePhenotype[plaquephenotype == -999] <- NA
AEDB$OverallPlaquePhenotype[plaquephenotype == -999] <- NA
AEDB$OverallPlaquePhenotype[plaquephenotype == 1] <- "fibrous"
AEDB$OverallPlaquePhenotype[plaquephenotype == 2] <- "fibroatheromatous"
AEDB$OverallPlaquePhenotype[plaquephenotype == 3] <- "atheromatous"
detach(AEDB)

# AEDB.temp <- subset(AEDB,  select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "plaquephenotype", "OverallPlaquePhenotype"))
# require(labelled)
# AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
# AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
# AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
# 
# DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
# 
# rm(AEDB.temp)

```

We will also fix and inverse-rank normal transform the continuous (manually)
scored plaque phenotypes.

```{r IRNT PlaquePhenotypes}
AEDB$macmean0 <- as.numeric(AEDB$macmean0)
AEDB$smcmean0 <- as.numeric(AEDB$smcmean0)
AEDB$neutrophils <- as.numeric(AEDB$neutrophils)
AEDB$Mast_cells_plaque <- as.numeric(AEDB$Mast_cells_plaque)
AEDB$vessel_density_averaged <- as.numeric(AEDB$vessel_density_averaged)

AEDB$MAC_rankNorm <- qnorm((rank(AEDB$macmean0, na.last = "keep") - 0.5) / sum(!is.na(AEDB$macmean0)))
AEDB$SMC_rankNorm <- qnorm((rank(AEDB$smcmean0, na.last = "keep") - 0.5) / sum(!is.na(AEDB$smcmean0)))
AEDB$Neutrophils_rankNorm <- qnorm((rank(AEDB$neutrophils, na.last = "keep") - 0.5) / sum(!is.na(AEDB$neutrophils)))
AEDB$MastCells_rankNorm <- qnorm((rank(AEDB$Mast_cells_plaque, na.last = "keep") - 0.5) / sum(!is.na(AEDB$Mast_cells_plaque)))
AEDB$VesselDensity_rankNorm <- qnorm((rank(AEDB$vessel_density_averaged, na.last = "keep") - 0.5) / sum(!is.na(AEDB$vessel_density_averaged)))

```

```{r IRNT PlaquePhenotypes: Visualisation}
library(labelled)
AEDB$Gender <- to_factor(AEDB$Gender)
library(patchwork)

p1 <- ggpubr::gghistogram(AEDB, "macmean0", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "% of macrophages (CD68)",
                    xlab = "% per region of interest", 
                    ggtheme = theme_minimal())

p2 <- ggpubr::gghistogram(AEDB, "MAC_rankNorm", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "% of macrophages (CD68)",
                   xlab = "% per region of interest\ninverse-rank normalized number", 
                    ggtheme = theme_minimal())

p1 | p2 

p1 <- ggpubr::gghistogram(AEDB, "smcmean0", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "% of smooth muscle cells (SMA)",
                    xlab = "% per region of interest", 
                    ggtheme = theme_minimal())

p2 <- ggpubr::gghistogram(AEDB, "SMC_rankNorm", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "% of smooth muscle cells (SMA)",
                   xlab = "% per region of interest\ninverse-rank normalized number", 
                    ggtheme = theme_minimal())

p1 | p2 


p1 <- ggpubr::gghistogram(AEDB, "neutrophils", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "number of neutrophils (CD66b)",
                    xlab = "counts per plaque", 
                    ggtheme = theme_minimal())

p2 <- ggpubr::gghistogram(AEDB, "Neutrophils_rankNorm", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "number of neutrophils (CD66b)",
                   xlab = "counts per plaque\ninverse-rank normalized number", 
                    ggtheme = theme_minimal())

p1 | p2 


p1 <- ggpubr::gghistogram(AEDB, "Mast_cells_plaque", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "number of mast cells",
                    xlab = "counts per plaque", 
                    ggtheme = theme_minimal())

p2 <- ggpubr::gghistogram(AEDB, "MastCells_rankNorm", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "number of mast cells",
                   xlab = "counts per plaque\ninverse-rank normalized number", 
                    ggtheme = theme_minimal())

p1 | p2 


p1 <- ggpubr::gghistogram(AEDB, "vessel_density_averaged", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "number of intraplaque neovessels",
                    xlab = "counts per 3-4 hotspots", 
                    ggtheme = theme_minimal())

p2 <- ggpubr::gghistogram(AEDB, "VesselDensity_rankNorm", 
                    # y = "..count..", 
                    color = "white",
                    fill = "Gender",
                    palette = c("#1290D9", "#DB003F"), 
                    add = "median", 
                    #add_density = TRUE,
                    rug = TRUE,
                    #add.params =  list(color = "black", linetype = 2), 
                    title = "number of intraplaque neovessels",
                   xlab = "counts per 3-4 hotspots\ninverse-rank normalized number", 
                    ggtheme = theme_minimal())

p1 | p2 

rm(p1, p2)
```

Here we calculate the *plaque instability/vulnerability* index

```{r Plaque Vulnerability}
# Plaque vulnerability
require(labelled)
AEDB$Macrophages.bin <- to_factor(AEDB$Macrophages.bin)
AEDB$SMC.bin <- to_factor(AEDB$SMC.bin)
AEDB$IPH.bin <- to_factor(AEDB$IPH.bin)
AEDB$Calc.bin <- to_factor(AEDB$Calc.bin)
AEDB$Collagen.bin <- to_factor(AEDB$Collagen.bin)
AEDB$Fat.bin_10 <- to_factor(AEDB$Fat.bin_10)
AEDB$Fat.bin_40 <- to_factor(AEDB$Fat.bin_40)

table(AEDB$Macrophages.bin)
table(AEDB$Fat.bin_10)
table(AEDB$Collagen.bin)
table(AEDB$SMC.bin)
table(AEDB$IPH.bin)

# SPSS code

# 
# *** syntax- Plaque vulnerability**.
# COMPUTE Macro_instab = -999.
# IF macrophages.bin=2 Macro_instab=1.
# IF macrophages.bin=1 Macro_instab=0.
# EXECUTE.
# 
# COMPUTE Fat10_instab = -999.
# IF Fat.bin_10=2 Fat10_instab=1.
# IF Fat.bin_10=1 Fat10_instab=0.
# EXECUTE.
# 
# COMPUTE coll_instab=-999.
# IF Collagen.bin=2 coll_instab=0.
# IF Collagen.bin=1 coll_instab=1.
# EXECUTE.
# 
# 
# COMPUTE SMC_instab=-999.
# IF SMC.bin=2 SMC_instab=0.
# IF SMC.bin=1 SMC_instab=1.
# EXECUTE.
# 
# COMPUTE IPH_instab=-999.
# IF IPH.bin=0 IPH_instab=0.
# IF IPH.bin=1 IPH_instab=1.
# EXECUTE.
# 
# COMPUTE Instability=Macro_instab + Fat10_instab +  coll_instab + SMC_instab + IPH_instab.
# EXECUTE.

# Fix plaquephenotypes
attach(AEDB)
# mac instability
AEDB[,"MAC_Instability"] <- NA
AEDB$MAC_Instability[Macrophages.bin == -999] <- NA
AEDB$MAC_Instability[Macrophages.bin == "no/minor"] <- 0
AEDB$MAC_Instability[Macrophages.bin == "moderate/heavy"] <- 1

# fat instability
AEDB[,"FAT10_Instability"] <- NA
AEDB$FAT10_Instability[Fat.bin_10 == -999] <- NA
AEDB$FAT10_Instability[Fat.bin_10 == " <10%"] <- 0
AEDB$FAT10_Instability[Fat.bin_10 == " >10%"] <- 1

# col instability 
AEDB[,"COL_Instability"] <- NA
AEDB$COL_Instability[Collagen.bin == -999] <- NA
AEDB$COL_Instability[Collagen.bin == "no/minor"] <- 1
AEDB$COL_Instability[Collagen.bin == "moderate/heavy"] <- 0

# smc instability
AEDB[,"SMC_Instability"] <- NA
AEDB$SMC_Instability[SMC.bin == -999] <- NA
AEDB$SMC_Instability[SMC.bin == "no/minor"] <- 1
AEDB$SMC_Instability[SMC.bin == "moderate/heavy"] <- 0

# iph instability
AEDB[,"IPH_Instability"] <- NA
AEDB$IPH_Instability[IPH.bin == -999] <- NA
AEDB$IPH_Instability[IPH.bin == "no"] <- 0
AEDB$IPH_Instability[IPH.bin == "yes"] <- 1

detach(AEDB)

table(AEDB$MAC_Instability, useNA = "ifany")
table(AEDB$FAT10_Instability, useNA = "ifany")
table(AEDB$COL_Instability, useNA = "ifany")
table(AEDB$SMC_Instability, useNA = "ifany")
table(AEDB$IPH_Instability, useNA = "ifany")

# creating vulnerability index
AEDB <- AEDB %>% mutate(Plaque_Vulnerability_Index = factor(rowSums(.[grep("_Instability", names(.))], na.rm = TRUE)),
                                )

table(AEDB$Plaque_Vulnerability_Index, useNA = "ifany")

# str(AEDB$Plaque_Vulnerability_Index)

```

## Prepare baseline summary

We are interested in the following variables at baseline in the whole cohort.

-   Age (years)
-   Female sex (N, %)
-   Artery type (N, %)
-   Hospital (N, %)
<!-- -   Hypertension (N, %) -->
<!-- -   SBP (mmHg) -->
<!-- -   DBP (mmHg) -->
<!-- -   Diabetes mellitus (N, %) -->
<!-- -   Total cholesterol levels (mg/dL) -->
<!-- -   LDL cholesterol levels (mg/dL) -->
<!-- -   HDL cholesterol levels (mg/dL) -->
<!-- -   Triglyceride levels (mg/dL) -->
<!-- -   Use of statins (N, %) -->
<!-- -   Use of antiplatelet drugs (N, %) -->
<!-- -   BMI (kg/m²) -->
<!-- -   Smoking status (N, %) -->
<!--     -   Never smokers -->
<!--     -   Ex-smokers -->
<!--     -   Current smokers -->
<!-- -   History of CAD (N, %) -->
<!-- -   History of PAD (N, %) -->
<!-- -   Clinical manifestations -->
<!--     -   Asymptomatic -->
<!--     -   Amaurosis fugax -->
<!--     -   TIA -->
<!--     -   Stroke -->
<!-- -   eGFR (mL/min/1.73 m²) -->
<!-- -   stenosis -->
<!-- -   year of surgery -->
<!-- -   plaque characteristics -->

```{r Baseline AEDB: preparation}
cat("====================================================================================================\n")
cat("SELECTION THE SHIZZLE\n")

### Artery levels
# AEdata$Artery_summary: 
#           value                                                                                   label
# NOT USE - 0 No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA
# USE - 1                                                                  carotid (left & right)
# USE - 2                                               femoral/iliac (left, right or both sides)
# NOT USE - 3                                               other carotid arteries (common, external)
# NOT USE - 4                                   carotid bypass and injury (left, right or both sides)
# NOT USE - 5                                                         aneurysmata (carotid & femoral)
# NOT USE - 6                                                                                   aorta
# NOT USE - 7                                            other arteries (renal, popliteal, vertebral)
# NOT USE - 8                        femoral bypass, angioseal and injury (left, right or both sides)

### AEdata$informedconsent
#           value                                                                                           label
# NOT USE - -999                                                                                         missing
# NOT USE - 0                                                                                        no, died
# USE - 1                                                                                             yes
# USE - 2                                                             yes, health treatment when possible
# USE - 3                                                                        yes, no health treatment
# USE - 4                                                yes, no health treatment, no commercial business
# NOT USE - 5                                                          yes, no tissue, no commerical business
# NOT USE - 6                      yes, no tissue, no questionnaires, no medical info, no commercial business
# USE - 7                             yes, no questionnaires, no health treatment, no commercial business
# USE - 8                                          yes, no questionnaires, health treatment when possible
# NOT USE - 9                  yes, no tissue, no questionnaires, no health treatment, no commerical business
# USE - 10                               yes, no health treatment, no medical info, no commercial business
# NOT USE - 11 yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business
# USE - 12                                                     yes, no questionnaires, no health treatment
# NOT USE - 13                                                             yes, no tissue, no health treatment
# NOT USE - 14                                                               yes, no tissue, no questionnaires
# NOT USE - 15                                                  yes, no tissue, health treatment when possible
# NOT USE - 16                                                                                  yes, no tissue
# USE - 17                                                                     yes, no commerical business
# USE - 18                                     yes, health treatment when possible, no commercial business
# USE - 19                                                    yes, no medical info, no commercial business
# USE - 20                                                                          yes, no questionnaires
# NOT USE - 21                         yes, no tissue, no questionnaires, no health treatment, no medical info
# NOT USE - 22                  yes, no tissue, no questionnaires, no health treatment, no commercial business
# USE - 23                                                                            yes, no medical info
# USE - 24                                                  yes, no questionnaires, no commercial business
# USE - 25                                    yes, no questionnaires, no health treatment, no medical info
# USE - 26                  yes, no questionnaires, health treatment when possible, no commercial business
# USE - 27                                                      yes,  no health treatment, no medical info
# NOT USE - 28                                                                             no, doesn't want to
# NOT USE - 29                                                                              no, unable to sign
# NOT USE - 30                                                                                 no, no reaction
# NOT USE - 31                                                                                        no, lost
# NOT USE - 32                                                                                     no, too old
# NOT USE - 34                                            yes, no medical info, health treatment when possible
# NOT USE - 35                                             no (never asked for IC because there was no tissue)
# USE - 36                    yes, no medical info, no commercial business, health treatment when possible
# NOT USE - 37                                                                                    no, endpoint
# USE - 38                                                         wil niets invullen, wel alles gebruiken
# USE - 39                                           second informed concents: yes, no commercial business
# NOT USE - 40                                                                              nooit geincludeerd

cat("- sanity checking PRIOR to selection")
library(data.table)
require(labelled)
ae.gender <- to_factor(AEDB$Gender)
ae.hospital <- to_factor(AEDB$Hospital)
table(ae.gender, ae.hospital, dnn = c("Sex", "Hospital"))
ae.artery <- to_factor(AEDB$Artery_summary)
table(ae.artery, ae.gender, dnn = c("Sex", "Artery"))

rm(ae.gender, ae.hospital, ae.artery)

# I change numeric and factors manually because, well, I wouldn't know how to fix it otherwise
# to have this 'tibble' work with 'tableone'... :-)

AEDB$Age <- as.numeric(AEDB$Age)
AEDB$diastoli <- as.numeric(AEDB$diastoli)
AEDB$systolic <- as.numeric(AEDB$systolic)

AEDB$TC_finalCU <- as.numeric(AEDB$TC_finalCU)
AEDB$LDL_finalCU <- as.numeric(AEDB$LDL_finalCU)
AEDB$HDL_finalCU <- as.numeric(AEDB$HDL_finalCU)
AEDB$TG_finalCU <- as.numeric(AEDB$TG_finalCU)

AEDB$TC_final <- as.numeric(AEDB$TC_final)
AEDB$LDL_final <- as.numeric(AEDB$LDL_final)
AEDB$HDL_final <- as.numeric(AEDB$HDL_final)
AEDB$TG_final <- as.numeric(AEDB$TG_final)

AEDB$Age <- as.numeric(AEDB$Age)
AEDB$GFR_MDRD <- as.numeric(AEDB$GFR_MDRD)
AEDB$BMI <- as.numeric(AEDB$BMI)
AEDB$eCigarettes <- as.numeric(AEDB$eCigarettes)
AEDB$ePackYearsSmoking <- as.numeric(AEDB$ePackYearsSmoking)
AEDB$EP_composite_time <- as.numeric(AEDB$EP_composite_time)
AEDB$EP_major_time <- as.numeric(AEDB$EP_major_time)

require(labelled)
AEDB$Artery_summary <- to_factor(AEDB$Artery_summary)
AEDB$ORyear <- to_factor(AEDB$ORyear)
AEDB$Gender <- to_factor(AEDB$Gender)
AEDB$Hospital <- to_factor(AEDB$Hospital)
AEDB$KDOQI <- to_factor(AEDB$KDOQI)
AEDB$BMI_WHO <- to_factor(AEDB$BMI_WHO)
AEDB$DiabetesStatus <- to_factor(AEDB$DiabetesStatus)
AEDB$SmokerStatus <- to_factor(AEDB$SmokerStatus)
AEDB$AlcoholUse <- to_factor(AEDB$AlcoholUse)

AEDB$Hypertension.selfreport <- to_factor(AEDB$Hypertension1)
AEDB$Hypertension.selfreportdrug <- to_factor(AEDB$Hypertension2)
AEDB$Hypertension.composite <- to_factor(AEDB$Hypertension.composite)
AEDB$Hypertension.drugs <- to_factor(AEDB$Hypertension.drugs)

AEDB$Med.anticoagulants <- to_factor(AEDB$Med.anticoagulants)
AEDB$Med.all.antiplatelet <- to_factor(AEDB$Med.all.antiplatelet)
AEDB$Med.Statin.LLD <- to_factor(AEDB$Med.Statin.LLD)

AEDB$Stroke_Dx <- to_factor(AEDB$Stroke_Dx)
AEDB$CAD_history <- to_factor(AEDB$CAD_history)
AEDB$PAOD <- to_factor(AEDB$PAOD)
AEDB$Peripheral.interv <- to_factor(AEDB$Peripheral.interv)

AEDB$sympt <- to_factor(AEDB$sympt)
AEDB$Symptoms.3g <- to_factor(AEDB$Symptoms.3g)
AEDB$Symptoms.4g <- to_factor(AEDB$Symptoms.4g)
AEDB$Symptoms.5G <- to_factor(AEDB$Symptoms.5G)
AEDB$AsymptSympt <- to_factor(AEDB$AsymptSympt)
AEDB$AsymptSympt2G <- to_factor(AEDB$AsymptSympt2G)
AEDB$Symptoms.Update2G <- to_factor(AEDB$Symptoms.Update2G)
AEDB$Symptoms.Update3G <- to_factor(AEDB$Symptoms.Update3G)

AEDB$restenos <- to_factor(AEDB$restenos)
AEDB$stenose <- to_factor(AEDB$stenose)
AEDB$EP_composite <- to_factor(AEDB$EP_composite)
AEDB$EP_major <- to_factor(AEDB$EP_major)
AEDB$Macrophages.bin <- to_factor(AEDB$Macrophages.bin)
AEDB$SMC.bin <- to_factor(AEDB$SMC.bin)
AEDB$IPH.bin <- to_factor(AEDB$IPH.bin)
AEDB$Calc.bin <- to_factor(AEDB$Calc.bin)
AEDB$Collagen.bin <- to_factor(AEDB$Collagen.bin)
AEDB$Fat.bin_10 <- to_factor(AEDB$Fat.bin_10)
AEDB$Fat.bin_40 <- to_factor(AEDB$Fat.bin_40)
AEDB$OverallPlaquePhenotype <- to_factor(AEDB$OverallPlaquePhenotype)
AEDB$Plaque_Vulnerability_Index <- to_factor(AEDB$Plaque_Vulnerability_Index)

AEDB$Artery_summary <- to_factor(AEDB$Artery_summary)

AEDB$informedconsent <- to_factor(AEDB$informedconsent)

AEDB.full <- subset(AEDB,
                    informedconsent != "missing" & # we are really strict in selecting based on 'informed consent'!
                       informedconsent != "no, died" & 
                       informedconsent != "yes, no tissue, no commerical business" &
                       informedconsent != "yes, no tissue, no questionnaires, no medical info, no commercial business" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commerical business" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business" &
                       informedconsent != "yes, no tissue, no health treatment" &
                       informedconsent != "yes, no tissue, no questionnaires" &
                       informedconsent != "yes, no tissue, health treatment when possible" &
                       informedconsent != "yes, no tissue" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commercial business" &
                       informedconsent != "no, doesn't want to" &
                       informedconsent != "no, unable to sign" & 
                       informedconsent != "no, no reaction" & 
                       informedconsent != "no, lost" & 
                       informedconsent != "no, too old" & 
                       informedconsent != "yes, no medical info, health treatment when possible" &
                       informedconsent != "no (never asked for IC because there was no tissue)" &
                       informedconsent != "no, endpoint" & 
                       informedconsent != "nooit geincludeerd")
# AEDB.CEA[1:10, 1:10]
dim(AEDB.full)

```

```{r Baseline AEDB: creation}
cat("===========================================================================================\n")
cat("CREATE BASELINE TABLE\n")

# Baseline table variables
basetable_vars = c("Hospital", 
                   "Artery_summary",
                   "Age", "Gender",
                   # "ORyear", 
                   # "TC_finalCU", "LDL_finalCU", "HDL_finalCU", "TG_finalCU", 
                   # "TC_final", "LDL_final", "HDL_final", "TG_final", 
                   # "hsCRP_plasma",
                   # "systolic", "diastoli", "GFR_MDRD", "BMI", 
                   # "KDOQI", "BMI_WHO",
                   # "SmokerStatus", "AlcoholUse",
                   # "DiabetesStatus", 
                   # "Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs", 
                   # "Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD", 
                   # "Stroke_Dx", "sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
                   # "Symptoms.Update2G", 
                   # "Symptoms.Update3G",
                   # "restenos", "stenose",
                   # "CAD_history", "PAOD", "Peripheral.interv", 
                   # "EP_composite", "EP_composite_time", "EP_major", "EP_major_time",
                   # "MAC_rankNorm", "SMC_rankNorm", "Macrophages.bin", "SMC.bin",
                   # "Neutrophils_rankNorm", "MastCells_rankNorm",
                   # "IPH.bin", "VesselDensity_rankNorm",
                   # "Calc.bin", "Collagen.bin", 
                   # "Fat.bin_10", "Fat.bin_40", 
                   "OverallPlaquePhenotype" )
                   # , "Plaque_Vulnerability_Index")

basetable_bin = c("Hospital", 
                  "Artery_summary",
                  "Gender")
# basetable_bin

basetable_con = basetable_vars[!basetable_vars %in% basetable_bin]
# basetable_con
```

# Athero-Express Biobank Study Baseline Characteristics

Showing the baseline table of the whole Athero-Express Biobank.

```{r Baseline AEDB: Visualize AEDB}
# Create baseline tables
# http://rstudio-pubs-static.s3.amazonaws.com/13321_da314633db924dc78986a850813a50d5.html
AEDB.tableOne = print(CreateTableOne(vars = basetable_vars, 
                                         factorVars = basetable_bin,
                                         # strata = "Symptoms.4g",
                                         data = AEDB.full, includeNA = TRUE), 
                          nonnormal = c(), missing = TRUE,
                          quote = FALSE, noSpaces = FALSE, showAllLevels = TRUE, explain = TRUE, 
                          format = "pf", 
                          contDigits = 3)[,1:3]
```


## Match samples used

Here we match the full database with the samples used in **`EntropyMasker`**. 

```{r}
# Remove duplicate rows of the dataframe
library(dplyr)
temp <- EM_samples %>% select(., "Studynumber") %>%
  distinct(.)

AEDB_EM <- merge(temp,
                 AEDB,
                 by.x = "Studynumber",
                 by.y = "STUDY_NUMBER", 
                 sort = FALSE)

AEDB_EM.full <- subset(AEDB_EM,
                    informedconsent != "missing" & # we are really strict in selecting based on 'informed consent'!
                       informedconsent != "no, died" & 
                       informedconsent != "yes, no tissue, no commerical business" &
                       informedconsent != "yes, no tissue, no questionnaires, no medical info, no commercial business" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commerical business" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business" &
                       informedconsent != "yes, no tissue, no health treatment" &
                       informedconsent != "yes, no tissue, no questionnaires" &
                       informedconsent != "yes, no tissue, health treatment when possible" &
                       informedconsent != "yes, no tissue" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info" &
                       informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commercial business" &
                       informedconsent != "no, doesn't want to" &
                       informedconsent != "no, unable to sign" & 
                       informedconsent != "no, no reaction" & 
                       informedconsent != "no, lost" & 
                       informedconsent != "no, too old" & 
                       informedconsent != "yes, no medical info, health treatment when possible" &
                       informedconsent != "no (never asked for IC because there was no tissue)" &
                       informedconsent != "no, endpoint" & 
                       informedconsent != "nooit geincludeerd")
# AEDB_EM.full[1:10, 1:10]

rm(temp)
```



```{r Baseline AEDB: Visualize AEDB EM FULL}
# Create baseline tables
# http://rstudio-pubs-static.s3.amazonaws.com/13321_da314633db924dc78986a850813a50d5.html
AEDB_EM.full.tableOne = print(CreateTableOne(vars = basetable_vars, 
                                         factorVars = basetable_bin,
                                         # strata = "Symptoms.4g",
                                         data = AEDB_EM.full, includeNA = TRUE), 
                          nonnormal = c(), missing = TRUE,
                          quote = FALSE, noSpaces = FALSE, showAllLevels = TRUE, explain = TRUE, 
                          format = "pf", 
                          contDigits = 3)[,1:3]
```

```{r Baseline AEDB: Visualize AEDB EM}

basetable_vars = c("Hospital", 
                   "Artery_summary")

basetable_bin = c("Hospital", 
                  "Artery_summary")

# Create baseline tables
# http://rstudio-pubs-static.s3.amazonaws.com/13321_da314633db924dc78986a850813a50d5.html
AEDB_EM.tableOne = print(CreateTableOne(vars = basetable_vars, 
                                         factorVars = basetable_bin,
                                         # strata = "Symptoms.4g",
                                         data = AEDB_EM, includeNA = TRUE), 
                          nonnormal = c(), missing = TRUE,
                          quote = FALSE, noSpaces = FALSE, showAllLevels = TRUE, explain = TRUE, 
                          format = "pf", 
                          contDigits = 3)[,1:3]
```

## Baseline writing

Let's save the baseline characteristics of the Athero-Express Biobank Study.

```{r Baseline SampleSelection: write}
# Write basetable

require(openxlsx)

write.xlsx(as.data.frame(AEDB.tableOne), 
           file = paste0(BASELINE_loc, "/",Today,".",PROJECTNAME,".AE.BaselineTable.xlsx"), 
           rowNames = TRUE, 
           colNames = TRUE, 
           sheetName = "AE_Base", overwrite = TRUE)

write.xlsx(as.data.frame(AEDB_EM.tableOne), 
           file = paste0(BASELINE_loc, "/",Today,".",PROJECTNAME,".AE.EM.59.BaselineTable.xlsx"), 
           rowNames = TRUE, 
           colNames = TRUE, 
           sheetName = "AE_Base_EM_59", overwrite = TRUE)

write.xlsx(as.data.frame(AEDB_EM.full.tableOne), 
           file = paste0(BASELINE_loc, "/",Today,".",PROJECTNAME,".AE.EM.56.BaselineTable.xlsx"), 
           rowNames = TRUE, 
           colNames = TRUE, 
           sheetName = "AE_Base_EM_56", overwrite = TRUE)

```

We will also write the newly prepared AEDB selected for this study which we can use in downstream analyses. 

```{r}
saveRDS(AEDB_EM.full, file = paste0(OUT_loc, "/",Today,".",PROJECTNAME,".AEDB.EM.FULL.RDS"))
saveRDS(AEDB.full, file = paste0(OUT_loc, "/",Today,".",PROJECTNAME,".AEDB.FULL.RDS"))
saveRDS(AEDB, file = paste0(OUT_loc, "/",Today,".",PROJECTNAME,".AEDB.raw.RDS"))

```


# Session information

------------------------------------------------------------------------------------------------------------------------

    Version:      v1.0.3
    Last update:  2022-08-23
    Written by:   Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
    Description:  Script to get some Athero-Express Biobank Study baseline characteristics.
    Minimum requirements: R version 3.4.3 (2017-06-30) -- 'Single Candle', Mac OS X El Capitan

    **MoSCoW To-Do List**
    The things we Must, Should, Could, and Would have given the time we have.
    _M_

    _S_

    _C_

    _W_

    **Changes log**
    * v1.0.3 Updated baseline characteristics.
    * v1.0.2 Simplified the initial script. It now outputs the relevant R-objects (as .RDS).
    * v1.0.1 Update to main AEDB (there is an error in the Age-variable in the new version).
    * v1.0.0 Initial version. Add 'plaque vulnerability index', Fixed baseline table, added codes, and results. Major update to WORCS system.

------------------------------------------------------------------------------------------------------------------------

```{r eval = TRUE}
sessionInfo()
```

# Saving environment

```{r Saving}
save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".AEDB.EM.baseline.RData"))
```

+-----------------------------------------------------------------------------------------------------------------------------------------+
| <sup>© 1979-2022 Sander W. van der Laan | s.w.vanderlaan[at]gmail.com | [swvanderlaan.github.io](https://swvanderlaan.github.io).</sup> |
+-----------------------------------------------------------------------------------------------------------------------------------------+