correlation_analysis.Rmd

---
title: "Correlation 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{correlation_analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

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

# Load libraries
library(bnaiaR)
library(bnlearn)
library(tidyr)
library(dplyr)
library(ggplot2)
library(ggthemes)

# Save plots as files? 
# Warning: Bad layout with knitr. Works well if chunks are run without rendering.
SAVEPLOTS <- FALSE
PLOTPATH <- Sys.getenv("PLOTPATH")
PLOTFORMAT <- "png" # "svg" or "png"
PLOTWIDTH = 16
PLOTHEIGHT = 9
```

# Load Data

```{r}
dataraw <- exp11_dat$abndata
```

```{r}
dataraw %>%
  mutate(id = row_number()) %>%
  tidyr::pivot_longer(cols = Gender, names_to = "Gender", values_to = "count") %>%
  pivot_wider(names_from = c("Gender", "count"), values_from = "count", values_fn = function(x){
    if(as.numeric(x)){
      1
      } else if(is.na(x)){
        0}
    }
  )
```

# One-hot encode data

```{r}
split2 <- function(x= dataraw$Gender){
  vals <- unique(x)
  df <- as.data.frame(matrix(ncol = length(vals), dimnames = list(NULL, vals)))
  
  v <- 1
  while (v <= length(vals)) {
    for (i in 1:length(x)) {
      if (x[[i]] == vals[v]){
        df[i,v] <- 1 
      } else {
        df[i,v] <- 0
      }
    }
    v <- v+1
  }
  
  return(df)
}
split2()
```


```{r}
dfonehot <- lapply(dataraw, split2) %>%
  do.call("cbind",.)
str(dfonehot)
```


```{r}
dfonehot <- dfonehot %>%
  rename(c("Sex = Female" = "Gender.Female",
           "Sex = Male" = "Gender.Male",
           "Fam. History = No" = "Positive.famillial.history.No",
           "Fam. History = Yes" = "Positive.famillial.history.Yes",
           "Age = F" = "AgeDiag.group.F",
           "Age = G" = "AgeDiag.group.G",
           "Age = D" = "AgeDiag.group.D",
           "Age = E" = "AgeDiag.group.E",
           "Age = C" = "AgeDiag.group.C",
           "Age = B" = "AgeDiag.group.B",
           "Age = A" = "AgeDiag.group.A",
           "Hypertension = Never" = "Hypertension.Never",
           "Hypertension = AnyType" = "Hypertension.AnyType",
           "Smoking = Current" = "Smoking_Current_Former_No.Current",
           "Smoking = Former" = "Smoking_Current_Former_No.Former",
           "Smoking = No" = "Smoking_Current_Former_No.No",
           "Location = High risk" = "location.grouped.High",
           "Location = Medium risk" = "location.grouped.Medium",
           "Location = Low risk" = "location.grouped.Low",
           "IA size = A" = "IAsize.groups.merged.A",
           "IA size = B" = "IAsize.groups.merged.B",
           "IA size = C" = "IAsize.groups.merged.C",
           "Multiple IAs = Yes" = "Multiple.IAs.Yes",
           "Multiple IAs = No" = "Multiple.IAs.No",
           "Ruptured IA = Yes" = "Ruptured_IA.Yes",
           "Ruptured IA = No" = "Ruptured_IA.No")) %>% 
  select(c("Sex = Female", 
           "Sex = Male",
           "Fam. History = Yes",
           "Fam. History = No",
           "Age = A",
           "Age = B",
           "Age = C",
           "Age = D",
           "Age = E",
           "Age = F",
           "Age = G",
           "Hypertension = AnyType",
           "Hypertension = Never",
           "Smoking = Current",
           "Smoking = Former",
           "Smoking = No",
           "Location = High risk",
           "Location = Medium risk",
           "Location = Low risk",
           "IA size = A",
           "IA size = B",
           "IA size = C",
           "Multiple IAs = Yes",
           "Multiple IAs = No",
           "Ruptured IA = Yes",
           "Ruptured IA = No"))
```

# Chi-square Test


make empty adjacency matrix

chi square visualization:
for each variable pair:
  calculate chi-square
  if chi-square p-value < 0.05:
    variables are probably related
    write 1 in adjacency matrix
  else
    variables are probably unrelated
    write 0 in adjacency matrix
    
```{r}
df <- dataraw %>%
  rename(c("Sex" = "Gender",
           # "Sex = Female" = "Gender.Female",
           # "Sex = Male" = "Gender.Male",
           "Fam. History" = "Positive.famillial.history",
           # "Fam. History = No" = "Positive.famillial.history.No",
           # "Fam. History = Yes" = "Positive.famillial.history.Yes",
           "Age" = "AgeDiag.group",           
           # "Age = F" = "AgeDiag.group.F",
           # "Age = G" = "AgeDiag.group.G",
           # "Age = D" = "AgeDiag.group.D",
           # "Age = E" = "AgeDiag.group.E",
           # "Age = C" = "AgeDiag.group.C",
           # "Age = B" = "AgeDiag.group.B",
           # "Age = A" = "AgeDiag.group.A",
           "Hypertension" = "Hypertension",
           # "Hypertension = Never" = "Hypertension.Never",
           # "Hypertension = AnyType" = "Hypertension.AnyType",
           "Smoking" = "Smoking_Current_Former_No",
           # "Smoking = Current" = "Smoking_Current_Former_No.Current",
           # "Smoking = Former" = "Smoking_Current_Former_No.Former",
           # "Smoking = No" = "Smoking_Current_Former_No.No",
           "Location" = "location.grouped",
           # "Location = High risk" = "location.grouped.High",
           # "Location = Medium risk" = "location.grouped.Medium",
           # "Location = Low risk" = "location.grouped.Low",
           "IA size" = "IAsize.groups.merged",
           # "IA size = A" = "IAsize.groups.merged.A",
           # "IA size = B" = "IAsize.groups.merged.B",
           # "IA size = C" = "IAsize.groups.merged.C",
           "Multiple IAs" = "Multiple.IAs",
           # "Multiple IAs = Yes" = "Multiple.IAs.Yes",
           # "Multiple IAs = No" = "Multiple.IAs.No",
           "Ruptured IA" = "Ruptured_IA"
           # "Ruptured IA = Yes" = "Ruptured_IA.Yes",
           # "Ruptured IA = No" = "Ruptured_IA.No"
           ))
```

```{r}
varnames <- names(df)
varno <- length(varnames)
adjmat <- matrix(rep(0, varno^2), nrow = varno, dimnames = list(varnames, varnames))
adjmat

for(r in 1:varno){
  for(c in 1:varno){
    chi <-chisq.test(df[[varnames[r]]], df[[varnames[c]]])
    if (chi$p.value < 0.05){
          adjmat[r,c] <- 1
    } 
  }
}
adjmat
```


```{r}
# keep only lower diagonal
adjmat[upper.tri(adjmat)] <- NA
diag(adjmat) <- NA


plt <- adjmat %>%
  as.data.frame() %>%
  tibble::rownames_to_column("rowvars") %>%
  pivot_longer(-c(rowvars), names_to = "colvars", values_to = "chisq") %>%
  mutate(rowvars = forcats::fct_inorder(rowvars),
         colvars = forcats::fct_inorder(colvars)) %>%

  
  # ggplot(aes(x=forcats::fct_rev(colvars), y=rowvars, fill=chisq)) + 
  ggplot() +
  aes(x=colvars, y=forcats::fct_rev(rowvars), fill=as.factor(chisq))+
  # ggplot(aes(x=colvars, y=rowvars, fill=chisq)) + 
  geom_raster() +
  # scale_fill_continuous(na.value = "transparent") +
  # scale_fill_grey(na.value = "transparent")+
  # scale_fill_steps(low = "lightgray", high="black", na.value = "transparent",
                   # name = "Correlation")+
  # scale_fill_stepsn(colors = c("lightgray", "black"), values = c(0,1), na.value = "transparent")+
  scale_fill_manual(values = c("lightgray", "black"), na.value = "transparent",
                   # name = "Chi^2 p-value < 0.05", labels = c("No: Independent", "Yes: Dependent", " ")
                   name = expression(p-value < 0.05), labels = expression(H[0]~not~rejected, H[0]~rejected))+
  ggtitle(expression(Chi^{2}~test~of~independence))+
  ggthemes::theme_tufte()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title = element_blank(),
        legend.position = c(0.7,0.7),
        legend.text.align = 0)

if (SAVEPLOTS) {
  plt <- plt + 
    theme(axis.text = element_text(size = 30),
          legend.text = element_text(size = 30),
          legend.title = element_text(size = 32),
          title = element_text(size = 32))
  ggsave(
    plot = plt,
    filename = "chisquareplt.png",
    path = PLOTPATH,
    width = PLOTWIDTH, height = PLOTHEIGHT, dpi = 600
  )
  } else {
  plt
}
```


# Spearman's Correlation Matrix

```{r}
# mat : is a matrix of data
# ... : further arguments to pass to the native R cor.test function
cor.mtest <- function(mat, ...) {
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat<- matrix(NA, n, n)
    diag(p.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], method = "spearman", ...)
            p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
        }
    }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

dfcor <- cor(dfonehot, method = "spearman")
head(dfcor)

# matrix of the p-value of the correlation
p.mat <- cor.mtest(dfonehot)
head(p.mat[, 1:5])

# Specialized the insignificant value according to the significant level
corrplot::corrplot(dfcor, type="upper", order="original", p.mat = p.mat, sig.level = 0.01)
```


```{r}
# Leave blank on no significant coefficient
plt <- corrplot::corrplot(dfcor, type="upper", 
                   
                   method = "color",
                   # method = "shade",
                   # method = "ellipse",
                   
                   # addCoef.col = "black", # Add coefficient of correlation
                   
                   tl.col="black", tl.srt=45, tl.cex = 0.6, #Text label color, rotation, size
                   
                   # Combine with significance
                   p.mat = p.mat,
                   
                   # sig.level = 0.01, insig = "blank", 
                   insig = 'label_sig', sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9, pch.col = 'grey20',
                   
                   # order = "hclust", addrect = 1, rect.col = "blue", rect.lwd = 4,
                   order="original",
                   
                   # hide correlation coefficient on the principal diagonal
                   diag=F
                   )

if (SAVEPLOTS){
  PLOTNAME <- "corrplot_spearman"
  dev.print(png, filename =paste0(PLOTPATH, "/", PLOTNAME, ".png"),
            res = 600, width = 6000)
  dev.off()
  } 
```

# Cramer's V

```{r}
varnames <- names(dfonehot)
varno <- length(varnames)
cramervmat <- matrix(rep(0, varno^2), nrow = varno, dimnames = list(varnames, varnames))
cramervmat

# for(r in 1:varno){
#   for(c in 1:varno){
#     phi <- DescTools::CramerV(dfonehot[[varnames[r]]], dfonehot[[varnames[c]]], conf.level = 0.95)
#     cramervmat[r,c] <- phi[[1]]
#   } 
# }

cramervmat <- DescTools::PairApply(dfonehot, DescTools::CramerV)

cramervmat

cramervmat.disc <- cramervmat
for(r in 1:nrow(cramervmat.disc)){
  for(c in 1:ncol(cramervmat.disc)){
    if(cramervmat.disc[r,c] <= 0.1){
      cramervmat.disc[r,c] <- 0.1
    } else if((cramervmat.disc[r,c] > 0.1) & (cramervmat.disc[r,c] <= 0.3)){
      cramervmat.disc[r,c] <- 0.3
    } else if((cramervmat.disc[r,c] > 0.3) & (cramervmat.disc[r,c] <= 0.5)){
      cramervmat.disc[r,c] <- 0.5
    } else if((cramervmat.disc[r,c] > 0.5) & (cramervmat.disc[r,c] <= 1)){
      cramervmat.disc[r,c] <- 1
    }
  } 
}
```

```{r}
# Leave blank on no significant coefficient
plt <- corrplot::corrplot(cramervmat, type="upper", 
                   
                   method = "color",
                   # method = "shade",
                   # method = "ellipse",
                   
                   # addCoef.col = "black", # Add coefficient of correlation
                   
                   tl.col="black", tl.srt=45, tl.cex = 0.6, #Text label color, rotation, size
                   
                   # Combine with significance
                   # p.mat = p.mat,
                   
                   # sig.level = 0.01, insig = "blank", 
                   # insig = 'label_sig', sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9, pch.col = 'grey20',
                   
                   # order = "hclust", addrect = 1, rect.col = "blue", rect.lwd = 4,
                   order="original",
                   
                   # hide correlation coefficient on the principal diagonal
                   diag=F,
                   
                   is.corr = F
                   )

if (SAVEPLOTS){
  PLOTNAME <- "cramerVplot"
  dev.print(png, filename =paste0(PLOTPATH, "/", PLOTNAME, ".png"),
            res = 600, width = 6000)
  dev.off()
  } 
```

Chi-square says that there is a significant relationship between variables, but it does not say just how significant and important this is. Cramer's V is a post-test to give this additional information.

Cramer's V varies between 0 and 1. Close to 0 it shows little association between variables. Close to 1, it indicates a strong association.

Where the table is 2 x 2, use Phi.

Cramer's V is named after the Swedish mathematician and statistician Harald Cramer.
(source: http://changingminds.org/explanations/research/analysis/cramers_v.htm)

# Phi coefficient

```{r}
varnames <- names(dfonehot)
varno <- length(varnames)
phimat <- matrix(rep(0, varno^2), nrow = varno, dimnames = list(varnames, varnames))
phimat

for(r in 1:varno){
  for(c in 1:varno){
    phi <- DescTools::Phi(dfonehot[[varnames[r]]], dfonehot[[varnames[c]]])
    phimat[r,c] <- phi[[1]]
  } 
}

phimat <- DescTools::PairApply(dfonehot, DescTools::Phi)

phimat
```

```{r}
# Leave blank on no significant coefficient
plt <- corrplot::corrplot(phimat, type="upper", 
                   
                   method = "color",
                   # method = "shade",
                   # method = "ellipse",
                   
                   # addCoef.col = "black", # Add coefficient of correlation
                   
                   tl.col="black", tl.srt=45, tl.cex = 0.6, #Text label color, rotation, size
                   
                   # Combine with significance
                   # p.mat = p.mat,
                   
                   # sig.level = 0.01, insig = "blank", 
                   # insig = 'label_sig', sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9, pch.col = 'grey20',
                   
                   # order = "hclust", addrect = 1, rect.col = "blue", rect.lwd = 4,
                   order="original",
                   
                   # hide correlation coefficient on the principal diagonal
                   diag=F,
                   
                   is.corr = F
                   )

if (SAVEPLOTS){
  PLOTNAME <- "phiplot"
  dev.print(png, filename =paste0(PLOTPATH, "/", PLOTNAME, ".png"),
            res = 600, width = 6000)
  dev.off()
  } 
```